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A correlated  independent  particle  (CIP)  model  is  formulated  where  the  effective 
Hamiltonian  is  composed  of  the  Fock  operator  and  a correlation  potential.  The  CIP 
model  can  be  regarded  as  a generalized  Hartree-Fock  model  that  includes  correla- 
tion, or  as  an  approximation  to  density  matrix  functional  theory,  where  the  density 
matrix  is  restricted  to  be  idempotent.  Within  the  model,  the  kinetic  energy  and  the 
exchange  energy  can  be  expressed  exactly,  leaving  the  correlation  energy  functional 
as  the  remaining  unknown.  Inspired  by  the  extended  Koopman’s  theorem,  the  condi- 
tion of  exactness  of  the  ionization  potentials  and  electron  affinities  as  the  eigenvalues 
of  the  effective  one-particle  operator  is  chosen  on  the  correlation  potential. 

If  no  truncation  is  introduced,  the  equation-of-motion  coupled-cluster  method 
for  ionization  potentials  and  electron  affinities  (IP-EOM-CC/EA-EOM-CC)  yields 
exact  ionization  potentials  and  electron  affinities.  An  energy-dependent  correlation 
potential  is  obtained  from  a partitioned  equation-of-motion  approach  which  repro- 
duces IP-EOM-CC  and  EA-EOM-CC  results  as  orbital  energies.  As  a consequence  of 
the  energy  dependence  of  the  correlation  potential,  a different  effective  Hamiltonian 
corresponds  to  each  orbital  and  each  eigenvalue.  The  equivalence  of  the  correlation 
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potential  to  the  energy-dependent  self-energy  used  in  the  propagator  method  with 
respect  to  the  eigenvalues  is  shown  if  the  coupled-cluster  wave  function  is  used  for 
the  ground  state. 

The  energy  dependence  is  regarded  as  a disadvantage,  but  can  be  removed  by 
using  the  Fock  space  coupled-cluster  method  in  the  one-hole  and  the  one-particle 
space  (which  is  equivalent  to  the  IP-EOM-CC  and  EA-EOM-CC  method)  to  extract 
a correlation  potential.  The  resultant  correlation  potential  is  energy-independent 
and  universal  for  all  orbitals.  In  a one-step  procedure,  the  exact  FS-CC  results  for 
ionization  potentials  and  electron  affinities  are  obtained  within  the  CIP  model. 

Further,  the  infinite-order  correlation  potential  derived  from  the  Fock  space 
coupled-cluster  approach  is  approximated  by  a potential  correct  through  second- 
order  in  perturbation.  To  improve  the  orbitals  compared  to  the  Hartree-Fock 
method,  the  Brillouin-Brueckner  condition  correct  through  second-order  in  pertur- 
bation is  superimposed.  Results  for  ionization  potentials,  and  electron  affinities  and 
properties  are  given  for  some  sample  molecules. 
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CHAPTER  1 
INTRODUCTION 

Because  of  their  conceptional  simplicity,  independent  particle  models  play  an 
essential  role  in  electronic  structure  theory.  Their  solutions  (molecular  orbitals) 
provide  the  conceptional  framework  for  much  of  chemistry.  Different  independent 
particle  models  emphasize  different  criteria  in  their  generation  and  their  solutions; 
and  associated  eigenvalues  form  the  basis  for  elements  ranging  from  photoelectronic 
spectra  to  frontier  molecular-orbital  theories. 

The  most  widely  used  independent  particle  models  are  density  functional  the- 
ory and  Hartree-Fock  theory.  Theoretical  justification  for  the  former  was  given  by 
Hohenberg  and  Kohn  [1].  It  is  a density-based  theory  that  is  formally  exact  and 
incorporates  correlation  effects.  Today’s  popular  version  of  DFT  was  initialized  by 
Kohn  and  Sham  [2] . In  it  orbitals  were  introduced  to  describe  the  kinetic  energy  and 
the  exact  density.  In  the  resultant  model,  eigenvalue  equations  for  one-electron  func- 
tions and  their  corresponding  energies  are  solved.  Correlation  is  included  through 
an  unknown  (but  often  approximated  [3,  4])  potential.  This  potential  also  has  to  ac- 
count for  the  kinetic  energy  correction,  because  the  kinetic  energy  in  the  Kohn  Sham 
model  is  that  of  a noninteracting  system.  Since  the  functional  dependence  of  the  ex- 
change energy  on  the  density  is  not  known  either,  exchange  potentials  can  also  only 
be  approximated  [5-7].  The  main  disadvantage  in  DFT  is  that  the  exchange  corre- 
lation potential  has  to  describe  different  kinds  of  energy  corrections  that  even  scale 
differently  [8].  That  makes  the  functional  and  the  potential  difficult  to  approximate, 
and  systematic  improvements  of  the  various  potentials  are  generally  not  available. 
Recently,  exact  local  exchange  potentials  and  second-order  ab  initio  perturbation 
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correlation  potentials  were  developed  [9-11]  that  are  orbital-dependent  and  guar- 
antee convergence  to  the  exact  result  in  the  correlation  and  basis  set  limit.  Some 
of  the  traditional  problems  in  DFT  (like  the  self  interaction,  the  wrong  asymptotic 
behavior  of  the  potential,  the  integer-discontinuity  condition,  and  even  the  inability 
to  obtain  weak  interactions)  are  removed  by  the  use  of  orbital-dependent  potentials 
in  ab  initio  dft. 

The  second  well-known  independent  particle  model  is  the  Hartree-Fock  theory. 
It  is  a wave-function  theory  that  uses  the  variational  principle  to  determine  the 
energetically  best  wave-function  given  by  a single  determinant  [12].  From  the  vari- 
ation, a set  of  eigenvalue  equations  for  the  orbitals  and  the  corresponding  orbital 
energies  is  obtained.  By  making  the  ansatz  of  a determinant  of  orbitals  for  the 
wave-function,  the  Pauli  principle  is  fulfilled  and,  therefore,  within  the  model,  the 
exchange  is  treated  exactly.  Similarly,  within  this  independent  particle  model,  the 
kinetic  energy  is  exact.  Koopmans’  theorem  [13]  also  gives  meaning  to  all  the  orbital 
energies,  while  Kohn  Sham  theory  only  assigns  the  highest  orbital  energy  to  an  ion- 
ization potential.  But  correlation  is  not  included  in  Hartree-Fock  theory.  Correlation 
corrections  can  be  added  to  the  wave-function  and  the  energy  through  many-body 
perturbation  theory  [14,  15],  coupled-cluster  theory  [16]  and  configuration  interac- 
tion [17]  for  instance,  at  the  cost  of  losing  the  simplicity  of  the  independent  particle 
wave-function. 

An  alternative  route  to  an  exact  one-electron  picture  is  based  on  electron- 
propagator  theory  [18,  19].  The  Dyson  equation  leads  to  a pseudo-eigenvalue  equa- 
tion for  energy-dependent  orbitals  (Dyson  orbitals).  The  corresponding  one-electron 
operator  has  an  energy-independent  Hartree-Fock  part,  and  an  energy-dependent 
part  (namely  the  self-energy  that  accounts  for  correlation  effects  [20]).  The  self- 
energy can  be  expanded  in  different  orders  of  perturbation  theory.  Therefore,  the 
theory  provides  a way  to  systematic  improvement.  Independent  particle  models 
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based  on  ab  initio  electron-propagator  methods  have  been  successfully  applied  to 
calculate  ionization  potentials  and  Dyson  orbitals  for  different  systems  [21-23].  To 
each  ionization  potential  there  corresponds  a Dyson  orbital  and  a self-energy,  de- 
pendent on  that  ionization  potential  (which  makes  it  necessary  to  solve  the  pseudo- 
eigenvalue problem  for  each  ionization  potential  and  each  orbital  separately,  making 
the  Dyson  orbitals  an  overdetermined,  non-orthogonal  set).  However,  their  diag- 
onalization  leads  to  natural  orbitals  [24],  The  same,  of  course,  holds  for  electron 
affinities  and  the  corresponding  orbitals. 

Gilbert’s  [25]  extensions  to  the  Hohenberg-Kohn  theorems  made  the  recent  for- 
mulation of  density  matrix  functional  theory  possible.  Natural  orbital  function- 
als were  developed  in  which  the  density  matrix  is  expressed  in  natural  orbitals 
and  corresponding  occupation  numbers.  Natural  orbital  functionals  can  be  di- 
vided into  corrected  Hartree  models  and  corrected  Hartree-Fock  [26]  models.  Buijse 
and  Baehrends  [27]  developed  a corrected  Hartree  scheme,  which  was  extended  by 
Godecker  and  Umrigar  [28]  to  include  an  additional  correction  term  for  the  self 
interaction.  Though  extremely  simple,  natural  orbital  functionals  show  surprising 
success  in  calculating  correlation  energies  for  atoms  and  small  molecules.  Analyzed 
by  various  authors  [29-33]  the  major  drawbacks  of  such  models  are  overcorrelation, 
inability  to  describe  the  limiting  case  of  the  high-density  homogeneous  electron  gas, 
and  violation  of  the  N-representability  condition  of  the  implicitly  defined  second- 
order  reduced  density  matrix.  Other  authors  proposed  density  matrix  functionals 
based  on  the  contracted  Schroedinger  equation  [34-37]  and  the  hypervirial  theorem 
[38,  39]. 

A highly  desirable  model  would  be  a generalized  independent  particle  model 
that  retains  the  advantage  of  the  Hartree-Fock  theory  in  being  able  to  describe  the 
kinetic  energy  and  the  exchange  energy  exactly,  but  which  incorporates  correla- 
tion effects.  This  could  be  achieved  by  applying  the  Brillouin-Brueckner  condition 
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or  equivalently,  the  maximum  overlap  condition  that  requires  that  a single  deter- 
minant of  Brueckner  orbitals  has  maximum  overlap  with  the  exact  wave-function. 
The  Brueckner  method  originated  in  nuclear  physics  [40]  and  became  known  in 
some  quantum  chemistry  circles  from  its  generalization  by  Lowdin  to  an  exact  self 
consistent  field  theory  [41].  Larsson  [42,  43]  developed  a Hartree-Fock-like  model 
that  includes  correlation  using  the  Brillouin- Brueckner  condition.  Chiles  and  Dyk- 
stra  [44],  Purvis  and  Bartlett  [45]  and  Adamowicz  and  Bartlett  [46],  and  Handy 
et  al.  [47]  obtained  Brueckner  orbitals  in  coupled-cluster  theory.  Stolarczyk  and 
Monkhorst  formulated  an  effective  Hamiltonian  that  has  the  form  of  the  Fock  op- 
erator plus  a correlation  potential  [48]  by  applying  the  maximum  overlap  condition 
to  a coupled-cluster  reference  wave-function.  Scuseria  [49,  50]  implemented  the 
Brueckner  coupled-cluster  method,  using  the  effective  Brueckner  Hamiltonian  in  a 
self  consistent  scheme.  Recently,  Lindgren  and  Salomonson  [51]  made  valuable  con- 
tributions to  the  field.  Minimizing  the  energy  with  respect  to  the  partitioning  of  the 
Hamiltonian  lead  them  to  a Brueckner  effective  Hamiltonian,  which  they  expressed 
in  different  orders  perturbation. 

Our  goal  was  to  formulate  a correlated  independent  particle  (CIP)  model  that 
has  the  general  form  (/  + vc)(f)p  = £p(j)p,  where  / is  either  the  Fock  operator  or 
a Fock-like  operator,  and  vc  is  the  correlation  potential  to  be  determined,  while 
the  (j)p  are  the  one-particle  solutions.  However,  the  question  is  what  condition  is 
reasonable  to  impose  to  obtain  vc.  We  chose  to  impose  the  condition  of  exactness  of 
the  ionization  potentials  (IP’s)  and  electron  affinities  (EA’s)  as  the  negative  of  the 
orbital  energies,  ev.  This,  of  course,  is  closely  related  to  Dyson  theory,  if  the  potential 
has  an  energy-dependent  form  as  is  the  case  for  potentials  constructed  in  Section 
3.1.3.  This  choice  is  justified  by  the  extended  Koopmans’  theorem  (EKT)  in  which 
the  eigenvalues  of  two  distinct  effective  one-particle  operators  represent  the  IP’s  and 
the  EA’s  respectively.  The  EKT  method  was  independently  developed  by  Day,  Smith 
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and  Garrod  [52]  for  single  ionization  and  attachment  processes  and  a generalization 
of  it  for  ionization  that  includes  double  and  higher  ionizations  by  Morrell,  Parr, 
and  Levy  [53].  Several  authors  argued  [53-55]  that  at  least  the  lowest  IP  is  exact 
if  the  exact  wave-function  is  used.  Numerical  support  was  given  by  Morrison  [56] 
and  Sundholm  and  Olsen  [57,  58].  The  energy-independent  effective  one-electron 
operator  formulated  in  Section  3.4  can  be  viewed  as  a realization  of  the  extended 
Koopmans’  theorem  Fock  operator  [59],  where  the  ground  state  wave-function  is 
represented  by  the  coupled-cluster  wave-function  and  where  the  assumption  that 
ionization  and  electron  capture  can  be  described  by  the  creation  or  annihilation  of 
a single  spin  orbital  is  relaxed,  and  additional  excitations  are  allowed. 

Our  study  explored  different  ways  to  obtain  a correlation  potential  for  a general- 
ized Hartree-Fock  model  that  includes  electron  correlation:  the  correlated  indepen- 
dent particle  (CIP)  model.  The  CIP  model  can  also  be  viewed  as  an  approximation 
to  density  matrix  functional  theory.  The  correlation  potential  is  chosen  such  that  the 
eigenvalues  of  the  effective  one-particle  method  represent  the  exact  IP’s  and  EA’s. 
Chapter  2 gives  basic  theorems  and  ideas  leading  to  the  formulation  of  the  CIP 
model.  The  Chapter  3 describes  different  realizations  of  the  correlation  potential 
yielding  exact  IP’s  and  EA’s  as  eigenvalues  of  an  effective  Hamiltonian.  Equiva- 
lence and  similarities  among  the  methods  used  to  obtain  the  correlation  potential 
(i.e.  the  equation-of-motion  coupled-cluster  method,  the  Fock  space  coupled-cluster 
method  and  the  coupled-cluster  Green’s  function  method)  are  shown.  In  Chapter  4 
the  energy-independent  correlation  potential  constructed  in  Section  3.4  is  approxi- 
mated through  second-order  in  perturbation,  and  results  are  given  for  some  sample 
molecules.  To  improve  the  orbitals  of  the  CIP  model  compared  to  Hartree-Fock  or- 
bitals, an  additional  condition  is  applied:  the  Brillouin-Brueckner  condition.  Results 
for  dipole  moments  of  sample  molecules  using  those  orbitals  are  listed  in  Chapter  4. 
Chapter  5 presents  the  conclusions. 


CHAPTER  2 

CORRELATED  INDEPENDENT  PARTICLE  MODEL  AS  AN 
APPROXIMATION  TO  DENSITY  MATRIX  FUNCTIONAL  THEORY 

In  this  chapter  the  basic  ideas  and  theorems  used  and  developed  to  formulate 
the  correlated  independent  particle  (CIP)  model  are  outlined.  The  CIP  model  is  an 
independent  particle  model  that  includes  correlation  through  a correlation  potential. 
The  eigenvalues  of  the  employed  effective  Hamiltonian  are  orbital  energies,  which 
are  associated  with  formally  exact  IP’s  and  EA’s.  The  CIP  model  is  viewed  as 
an  extension  to  the  Hartree-Fock  model  or  as  an  approximation  to  density  matrix 
functional  theory  where  the  first-order  density  matrix  is  restricted  to  be  idempotent. 

Since  the  Hartree-Fock  theory  is  a special  case  of  the  CIP  model  and  density 
matrix  functional  theory,  the  Hartree-Fock  model  is  briefly  summarized  in  Section 
2.1.  In  Section  2.2  the  necessary  theorems  and  proofs  to  establish  the  existence  of  a 
universal  energy  functional  of  the  first-order  reduced-density  matrix  are  reported  and 
the  validity  of  the  variational  principle  is  shown.  Hartree-Fock  theory  is  shown  to 
be  a special  case  of  density  matrix  functional  theory.  In  Section  2.3  the  CIP  model 
is  introduced.  The  condition  of  exactness  of  the  ionization  potentials  (IP’s)  and 
electron  affinities  (EA’s)  as  the  eigenvalues  of  the  one-particle  operator  employed  in 
the  CIP  model  is  used  to  determine  the  correlation  potential.  This  choice  is  inspired 
by  the  extended  Koopmans’  theorem,  which  is  described  in  Section  2.4.  Finally  in 
Section  2.5  wesummarize  the  different  density  matrix  functionals  developed  in  the 
past. 
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2.1  Hartree-Fock  Theory 

The  Hartree-Fock  approximation  [12]  is  one  of  the  most  widely  used  indepen- 
dent particle  theories  in  quantum  chemistry.  It  is  a wave-function  method  attempt- 
ing to  solve  the  many-electron  problem  by  finding  an  approximate  solution  to  the 
nonrelativistic  time-independent  electronic  Schroedinger  equation.  The  simplicity 
of  the  wave-function  allows  for  a direct  chemical  interpretation.  The  eigenfunctions 
of  the  Fock  operator  are  orbitals  with  corresponding  orbital  energies  that  play  a 
vital  role  in  explaining  and  determining  chemical  reaction  mechanisms.  However, 
the  lack  of  accuracy  of  the  Hartree-Fock  approximation  prevents  the  prediction  of 
reaction  energies  and  other  electronic  properties  within  a few  kcal/mole.  To  reach 
higher  accuracy,  correlated  methods  like  many-body  perturbation  theory  [14,  15], 
coupled-cluster  theory  [16]  and  configuration  interaction  [17]  have  to  be  used.  But 
since  the  Hartree-Fock  approximation  usually  accounts  for  most  of  the  total  energy, 
it  is  often  used  as  a starting  point  for  correlated  methods. 

At  first  glance  the  Hartree-Fock  approximation  does  not  use  the  density  matrix 
as  its  basic  variable,  but  rather  the  wave-function.  However,  the  particular  ansatz  for 
the  wave-function  leads  to  eigenvalue  equations  that  are,  as  seen  later,  the  equations 
of  a special  case  in  density  matrix  functional  theory:  namely,  when  the  correlation 
energy  is  zero  (i.e.,  the  Hartree-Fock  energy  can  be  written  in  terms  of  the  first-order 
density  matrix  only). 

The  nonrelativistic  time-independent  Schroedinger  equation  is  given  by 

i/'F  = FhF,  (2.1) 

where  H is  the  Hamiltonian  operator  for  a system  of  nuclei  and  electrons  described 
by  position  vectors.  The  distance  between  the  ith  electron  and  the  Ath  nucleus 
is  riA,  the  distance  between  the  ith  electron  and  the  jth  electron  is  r^.  In  the 
Born-Oppenheimer  approximation,  the  electrons  move  in  a field  of  fixed  nuclei.  The 


8 


Hamiltonian  (in  Hartree  atomic  units)  is  then 


(2.2) 


N is  the  number  of  electrons  and  M is  the  number  of  nuclei.  Z&  is  the  atomic 
number  of  nucleus  A.  The  first  term  on  the  right-hand  side  is  the  operator  for  the 


kinetic  energy  of  the  electrons,  the  second  term  is  the  nuclear-electron  attraction 
operator,  and  the  third  term  is  the  operator  for  electron-electron  repulsion. 

The  variational  principle  states  that  the  exact  energy  is  always  lower  than  or 
equal  to  the  expectation  value  of  the  Hamiltonian  with  respect  to  a trial  normalized 
wave-function  $ 


for  the  wave-function  can  be  used  and  the  energy  minimized  to  obtain  the  best 
wave-function  of  that  particular  form. 

The  simplest  antisymmetric  wave-function  of  an  iV-electron  system  that  is  com- 
posed of  one-electron  functions  is  a Slater  determinant 


The  one-electron  functions  fa  depend  on  the  space-spin  coordinates  Xi  of  the  elec- 
trons. For  orthonormal  spin  orbitals,  the  expectation  value  of  the  Hamiltonian  with 
respect  to  the  Slater  determinant  is  given  by 


£exact  < (*\H\*)- 


(2.3) 


The  Dirac  notation  [60]  is  used  to  denote  integral  expressions.  A reasonable  ansatz 


fa(xi)  & (xi)  ...  0jv(xi) 

_ . ,1  01  (X2)  02  (x2)  •••  0v(x2) 

$0(Xi,X2,...,Xtf)  = 


0l(Xjv)  02(X;v)  •••  0Af(Xjv) 


(2.4) 
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h is  the  one-electron  operator  part  of  the  Hamiltonian  and  includes  the  kinetic 
energy  operator  and  the  nuclear-electron  attraction.  The  double  bar  1 1 denotes  the 
sum  of  the  coulomb  and  the  exchange  contributions  to  the  two-electron  integral 

(<t>i4>j\\fa<t>j)  = (fafalfafa)  ~ {fafa\fafa)- 

In  the  Hartree-Fock  model,  the  spin  orbitals  fa  are  varied  under  the  orthogo- 
nality constraint  (fa\cf)j)  — Sij  to  obtain  the  Slater  determinant  that  corresponds  to 
the  lowest  energy.  The  problem  of  minimizing  a functional  under  constraints  can  be 
solved  by  the  Lagrangian  multiplier  method.  The  functional  to  be  minimized  is 

N 

= Eq  — ^ ] eji((0t|0j)  ~ &ij ) ) (2.5) 

ij 


where  the  are  the  Lagrange  multipliers.  Variation  of  the  orbitals  leads  to  the 
Hartree-Fock  equations 

N 


fifa)  = 

j 


(2.6) 


/ is  the  Fock  operator,  defined  as 

N 

fix i)  = h{xi)  + E / ^(*2)^(1  - V12)faix2)dx2.  (2.7) 

3 

The  operator  V\2  permutes  electron  one  and  electron  two. 

For  a single  determinant  wave- function,  any  expectation  value  is  invariant  under 
a unitary  transformation  of  the  spin  orbitals.  Since  e is  a Hermitian  matrix,  it  is  al- 
ways possible  to  find  a unitary  transformation  that  diagonalizes  e.  The  transformed 
Hartree-Fock  equations 


f\fa)  = 6i\fa)  (2.8) 

are  known  as  canonical  Hartree-Fock  equations.  The  eigenfunctions  of  the  canonical 
Hartree-Fock  equations  are  the  canonical  spin  orbitals,  and  the  eigenvalues  are  the 
orbital  energies. 
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The  Fock  operator  has  a functional  dependence  on  the  occupied  spin  orbitals 
but  once  the  spin  orbitals  are  known,  the  Fock  operator  becomes  a well-defined  Her- 
mitian  operator  with  an  infinite  number  of  eigenfunctions.  The  N orbitals  with  the 
lowest  orbital  energies  are  the  occupied  orbitals  in  <f>0  labeled  i,  j, . . the  remaining 
spin  orbitals  are  the  virtual  orbitals  with  indices  a,  b, . . ..  Labels  p,q,...  refer  to 
general  orbitals.  The  orbital  energy  is  given  by 

£p  = (^pl/l^p) 

N 

= <^|/#P>  + £<^||^>,  (2-9) 

i 

which,  according  to  Koopmans’  theorem  [13],  is  equal  to  the  frozen  orbital  estimate 
of  the  negative  of  the  ionization  potential  (IP)  for  occupied  orbitals  and  similar  for 
the  negative  of  the  electron  affinity  (EA)  for  unoccupied  orbitals 

IP  = N~lEk  - nE0  = —ek  (2.10) 

EA  = nE0  - n+1Ec  = -ec.  (2.11) 

N~1Ek  is  the  expectation  value  of  the  Hamiltonian  with  respect  to  <f>0  where  an 
electron  in  orbital  k is  removed.  N+l Ec  is  the  expectation  value  of  the  Hamiltonian 
with  respect  to  <f>0  where  an  electron  in  orbital  c is  added. 

Adding  the  orbital  energies  Cj  of  Eq.  2.9  for  all  the  occupied  orbitals,  gives 

N N N 

Ylei  = + (2.i2) 

* * ij 

Comparing  Eq.  2.4  with  Eq.  2.12,  it  is  apparent  that  the  total  Hartree-Fock  energy 
is  not  the  sum  of  the  orbital  energies.  The  sum  of  orbital  energies  counts  the 
electron-electron  interaction  twice,  and  has  to  be  corrected  by  a factor  of  | to  yield 
the  correct  total  energy. 

For  molecules  with  an  even  number  of  electrons,  it  is  found  that  the  spin  orbitals 
are  essentially  pairwise  degenerate.  Restricted  Hartree-Fock  theory  takes  that  fact 
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explicitly  into  account.  A restricted  set  of  spin  orbitals  has  the  form 


{<pp(r)a(u) 
ipp(r)p(u) 


(2.13) 


where  the  (pp(r)  are  spatial  orbitals  dependent  on  the  spatial  coordinates  r.  a(ui) 
and  /3(lj)  are  spin  eigenfunctions  dependent  on  the  spin  coordinate  u,  fulfilling 
= (I3\P)  = 1 and  (a\P)  = (P\a)  = 0. 

To  convert  the  Hartree-Fock  equations  into  restricted  Hartree-Fock  equations, 
the  spin  integration  is  carried  out,  leading  to 


/ Wi)  = f.lv.), 


(2.14) 


with 

N 

f(r i)  = h(n)  + f ^(r2)rf21(2  - Vi2)ipj(r2)dr2.  (2.15) 

j 

The  closed-shell  orbital  energies  are 

N 

2 

= (<Pp\h\<Pp)  T ^ ^(2  ((flp(fii\ipp(pi)  — (</?p(^i|(pj(pp)),  (2.16) 

i 

for  which  Koopmans’  theorem  applies  as  well.  The  total  energy  is  given  by 

N_  N 

2 2 

Eo  = 2^2(ipi\h\ipi)  + ^(2  (<pj<pi\<pjipi)  - (<Pj<pi\<pi<pj)).  (2.17) 

i ij 


2.2  Hohenberg-Kohn-Gilbert  Theorem 

The  Hohenberg-Kohn-Gilbert  theorem  [25]  is  the  foundation  of  recent  devel- 
opments in  density  matrix  functional  theory.  Hohenberg  and  Kohn  [1]  proved  the 
existence  of  a universal  energy  functional  of  the  external  potential  and  the  particle 
density  that  is  a minimum  for  the  true  particle  density.  The  proof  was  given  for  a 
local  external  potential  v(x,  x ')  = 5(x  — x')v(x).  Gilbert  [25]  extended  the  theorem 
to  a nonlocal  external  potential.  He  showed  that  the  universal  energy  functional  of 
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the  nonlocal  external  potential  depends  on  the  one-particle  density  matrix  instead 
of  the  density. 

2.2.1  Energy  Functional  of  the  Density  Matrix 

The  von  Neumann  N-particle  statistical  density  kernel  is  given  by 

DN(x i . ,.xN,x[  . . .x'N)  = '^2/wi^i{x i . . .x/v)T*(x'1 . . .x'N),  (2.18) 

i 

where  the  'It,  are  a complete  set  of  orthonormal  N-particle  wave-functions,  the  wl 
are  positive  numbers  normalized  to 


^2wi  = l.  (2.19) 

i 

The  ensemble  N-representable  one-particle  density  kernel  can  be  obtained  from  an 
N-particle  kernel  by  contraction 

»>(*,*')  = nJ  Dn(x,  X2  . . . Xjv;  x1 , X2  ■ ■ ■ xu)dx2  . . . dxN . (2.20) 


Similarly,  the  ensemble  N-representable  two-particle  density  kernel  can  be  defined 
D2(x1x2,  x[x’2)  = 


N ' 


Z / 


J Dn(xi,x2,  x3...  xN]x[ , x'2,x3  . ..xN)dx3  . . . dxN, 

(2.21) 


and  so  on. 

Partition  the  Hamiltonian  in  the  following  way 


H = + + (2-22) 

i i i<j 

t is  the  kinetic  energy  operator,  v is  the  external  potential,  and  u is  the  electron- 
electron  interaction.  The  energy  can  be  written  as  a functional  of  the  N-particle 
density  kernel  and  the  external  potential  (t  and  u will  always  be  the  same;  their 
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functional  dependence  does  not  need  to  be  indicated) 

E{Dn,  v } = ( DnH ) = { DH ) + {D1v)  + ( D2u ).  (2.23) 

The  external  potential  v(x , x')  does  not  have  to  be  local,  but  is  restricted  to  the 
class  of  potentials  for  which  E is  real  and  has  a lower  bound.  The  ground  state  \&0 
is  assumed  to  be  nondegenerate,  so  that  Dq  = ToTq. 

From  Eq.  2.23  it  can  be  seen  that  the  energy  is  trivially  a function  of  the 
two-matrix.  In  1951  Coleman  [61]  attempted  to  vary  the  energy  with  respect  to 
the  two-matrix  and  was  surprised  to  find  an  energy  for  the  lithium  atom  that  was 
20%  below  the  experimental  value.  The  reason  for  the  failure  of  the  variational 
principle  is  that  the  space  of  two-matrices  in  which  he  varied  was  too  large.  The 
condition  that  must  be  fulfilled  by  the  two-matrices  is  that  they  correspond  to  a 
fermion  wave-function,  and  that  they  can  be  obtained  by  integrating  the  N-particle 
density  matrix  over  N — 2 particles.  Ensuring  this  condition  is  known  as  the  N- 
representability  problem,  currently  thought  to  be  an  unsolved  problem  (at  least 
practically)  for  two-matrices.  Therefore,  it  is  not  particularly  useful  to  write  the 
energy  as  a functional  of  the  two-matrix. 

In  Eq.  2.23,  the  ensemble  energy  functional  E = (DN H)  has  been  used  rather 
than  the  pure- state  energy  functional  E — (T|i/|\k),  because  the  constraints  needed 
to  apply  the  variational  principle  are  known  only  for  ensembles.  Next  it  is  shown 
that  the  same  stationary  values  are  obtained  for  ensemble  or  pure-state  energy  func- 
tionals. 

To  derive  the  stationary  conditions  for  ensembles,  the  constraints  (dqlTj)  = 
( i\j ) = Sij  and  J2wi  — 1 must  be  imposed  by  the  method  of  Lagrangian  multipliers. 
The  Wi  are  set  to  be  cos29i  to  satisfy  the  constraint  0 < Wi  < 1.  The  functional  to 
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be  varied  with  respect  to  'Iq,  4/}  and  0*  is 

£ = e b’>  + A(i  - X^)’  (2.24) 

hj 

where  A and  \ZJ  are  Lagrangian  multipliers.  The  variation  leads  to  the  following 
equations 


se 

ss 

6Vi 

dS_ 

d9i 


wtH^i  - Tj  A tj  = 0 

3 

3 

sm26i(Ei  — A)  = 0, 


(2.25) 

(2.26) 
(2.27) 


where  Ei  = (i\H\i).  Eq.  2.25  and  Eq.  2.26  reduce  to  the  single  condition 


(wi  - Wj)(i\H\j)  = 0. 


(2.28) 


From  Eq.  2.27  it  can  be  concluded  that  one  of  the  following  conditions 


( i ) Wi  = 1 (ii)  Ei  = A (Hi)  Wi  = 0 

must  be  satisfied.  If  ( i ) is  true,  then  all  other  Wj  — 0,  and  since  (i\H\j)  = 0 for 
i ^ j,  according  to  Eq.  2.28.  Therefore, 

H^i  = E^!i.  (2.30) 

If  (ii)  is  true,  then  all  states  for  which  Wi  ^ 0 must  belong  to  the  same  degenerate 
state,  and  Eq.  2.30  is  still  satisfied.  In  all  other  cases,  wt  = 0.  This  shows  that 
stationary  values  for  the  ensemble  energy  functional  are  the  same  as  for  the  pure- 
state  energy  functional. 

Extending  the  Hohenberg-Kohn  theorem  to  nonlocal  potentials  [25]  requires 
only  slight  modifications  of  the  proof  given  by  Hohenberg  and  Kohn  [1],  v1  and  v2 
are  two  external  potentials  that  give  rise  to  to  two  distinct  nondegenerate  ground 
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states  and  The  corresponding  Schrodinger  equations  are  given  by 


H1^  = Efol  (2.31) 

H2^2  = Efol  (2.32) 

The  reduced  one-particle  density  kernels  for  the  two  different  wave-functions  are 

‘nj  = A W)  ‘A?  = *J*J'  (2.33) 

2Z>J  = jV(2D,f)  2Z)J'  = tf2^’.  (2.34) 

Applying  the  variational  principle  leads  to 

El  = {lD"Hl)  < E1  — {2D*Hl)  (2.35) 

El  = {2D*H2)  <E2  = CD^H2),  (2.36) 


where  Eq  and  El  are  the  ground-state  energies.  E1  and  E2  are  necessarily  greater 
since  to  calculate  El  the  ground-state  density  matrix  of  system  2 (which  is  different 
from  the  ground-state  density  matrix  of  system  1)  is  used.  An  energy  difference  can 
be  defined 

A E = (E2  - El)  + (E1  - El)  > 0,  (2.37) 

where  Eq.  2.35  and  Eq.  2.36  have  been  used.  Rearranging  Eq.  2.37  and  inserting 
the  definitions  Eq.  2.35  and  2.36  leads  to 

A E = (E2  - El)  + (E1  - El)  (2.38) 

= C Dq  (H2  - H1))  + (2D"(Hl  - H2)).  (2.39) 

The  difference  between  the  two  Hamiltonians  is  the  difference  between  the  two  ex- 
ternal potentials  Sv  = v2  — v1,  which  is  by  definition  nonzero.  Then 
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A E = CD^Sv)  - (2D£5v)  (2.40) 

= ~(6Dq6v)  > 0,  (2.41) 

where  SDq  = 2D\  — 1D\. 

For  the  inequality  Eq.  2.41  to  hold,  8Dq  cannot  be  zero.  Therefore  if  Tq  and 
Tq  are  distinct,  then  their  one-particle  density  kernels  1D\  and  2Dq  must  also  be 
distinct,  implying  that  there  is  a one-to-one  correspondence  between  Dq  and  D\. 
Hereby,  the  existence  of  an  energy  functional  of  the  one-particle  reduced-density 
kernel  and  the  external  potential  is  established 

E{Dl,v]  = E{D^{D'0),v}  = <DJ»)  + F{£>J},  (2.42) 

where  F{Dq}  is  a universal  functional  of  D\  alone. 

The  domain  of  F{Dq}  (as  defined  in  Eq.  2.42)  consists  only  of  those  density 
kernels  that  can  be  constructed  from  the  nondegenerate  ground-state  wave-function 
in  a local  or  nonlocal  external  potential.  To  extend  this  domain,  define  a set  {DN} 
that  includes  Dq  and  any  other  N-particle  density  kernel  that  is  needed  for  a one- 
to-one  mapping  between  {DN}  and  the  set  of  one-particle  reduced-density  kernels 
{D1}.  The  set  {DN}  is  not  uniquely  defined,  in  particular  because  a set  {DN}  can 
be  found  for  which  a many-to-one  mapping  onto  {D1}  exists.  However,  here  it  is 
only  important  that  one  set  can  be  found,  which  allows  for  the  one-to-one  mapping. 
Then 

E{D\v}  = E{Dn{D1},  v}  = (D'v)  + FiD1},  (2.43) 

where  D 1 is  now  any  density  kernel  associated  with  a ground  state  in  an  external 
potential.  This  allows  for  the  definition  of  the  variational  principle 


E{Dl,v}  < E{D10  + 5D\v}. 


(2.44) 
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The  existence  of  a 1-matrix  functional  theory  was  also  shown  by  other  authors. 
Berrondo  and  Goscinski  [62]  added  a nonlocal  external  potential  to  the  N-particle 
Hamiltonian  and  then  obtained  a variational  principle  involving  the  1-matrix  for 
a local  external  potential  by  eliminating  the  nonlocal  external  source.  Donnelly 
and  Parr  [63]  showed  that  the  existence  of  a universal  variational  functional  of  the 
1-matrix  is  implied  in  the  original  Hohenberg-Kohn  theorem  [1]. 

2.2.2  V-Representability  Problem 

A limitation  of  F{D1}  is  that  it  is  undefined  for  any  D1  that  is  not  v-representable. 
A ground-state  v-representable  D 1 is  one  that  is  associated  with  an  antisymmetric 
ground-state  wave-function  of  some  Hamiltonian  with  a local  or  nonlocal  poten- 
tial. The  v-representability  constraint  is  very  severe,  since  most  one-particle  den- 
sity matrices  are  not  v-representable  (for  instance,  no  idempotent  one-matrix  is 
v-representable).  It  was  shown  by  Levy  [64]  with  the  constrained  search  formalism 
that  v-representability  is  not  required.  The  N-representability  constraint,  where  the 
one-matrix  must  be  derivable  from  an  antisymmetric  wave-function,  is  sufficient. 
The  proof  is  presented  in  the  following. 

Define  the  universal  functional 

W{D1}  = min(TjDi|u|\I>Di).  (2-45) 

W{D1}  searches  all  antisymmetric  wave-functions  T Di  that  yield  the  fixed  trial 
D1,  where  Dl  does  not  need  to  be  v-representable.  W{D1}  delivers  the  minimal 
expectation  value.  For  W{D1}  to  be  a valid  universal  variational  functional  for  any 
N- representable  Dl  two  theorems  need  to  be  proven. 


Theorem  I (D1v)  + (Dh)  + W{D1}  > E0 

Theorem  II  {D\v)  + {D\t)  + W{Dq}  = E0 


(2.46) 

(2.47) 
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Call  T"]1/1  that  wave-function  that  satisfies  the  right-hand  side  of  Eq.  2.45 

W{D1}  = (tf$>|tf£y*)  (2.48) 

and 

W{Dl0}  = (2.49) 

To  prove  Theorem  I,  insert  Eq.  2.48  into  the  left-hand  side  of  Eq.  2.46 


{D1v)  + ( DH ) + W{D1}  = + v + u\V™?).  (2.50) 

But  by  the  variational  principle 

{y™?\t  + v + u\V™?)  >E0.  (2.51) 

Adding  the  last  two  equations  concludes  the  proof  of  Theorem  I. 

To  prove  Theorem  II,  again  use  the  variational  principle 

E0  < (^f\t  + v + u\^)  (2.52) 

Eq  = (\Ho|t  + v + ^l^o)-  (2.53) 

From  there  it  is  seen  that 

(Dlv)  + (Dlt)  + (*0|«l*„>  < (Dfr)  + (D'0t)  + <^M»Sr>.  (2.54) 

which  leads  to 

(tfoMTo)  < (2.55) 

But  the  definition  of  4'™"  dictates  that 

U0 

(tfoM^o)  > (^|«|^>.  (2.56) 

The  last  two  equations  can  only  hold  at  the  same  time  if 

('PoM'I'o)  = |«|*S}">  = W{Dl). 


(2.57) 
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The  ground-state  energy  is  given  by 


Eo  — <*o|f  + v + w|Tq)  — (D^t)  + (Dqv)  + (tf0|u|tfo>.  (2.58) 


Substitution  of  Eq.  2.57  into  Eq.  2.58  completes  the  proof  of  Theorem  II.  Notice 
that  FID1}  = W{D1}  if  Dl  is  v-representable. 

In  addition  to  the  v-representability  constraint,  the  condition  that  the  ground 
state  is  nondegenerate  is  lifted  within  the  constrained  search  formalism.  If  T0  is 
degenerate,  then  all  of  the  ground-state  wave-functions  may  be  obtained,  one  at  a 
time,  by  the  foregoing  procedure. 

2.2.3  Euler  Equations 

Euler  equations  that  determine  the  density  matrix  can  now  be  derived  by  apply- 
ing the  variational  principle.  According  to  Coleman  [65],  a Hermitian  one-electron 
operator  D1(x,x')  will  be  ensemble  N-representable  if  and  only  if  f D1(x,x)  = N 
and  all  eigenvalues  satisfy  the  inequality  0 < np  < 1, 


The  {cj)i}  are  the  natural  orbitals  [66],  the  eigenfunctions  of  the  one-particle  density 
kernel.  The  N-represent ability  conditions  on  the  pure-state  one-particle  density 
kernel  are  not  known.  However,  as  already  shown  above,  one  obtains  the  same 
stationary  conditions  for  unreduced-ensemble  or  pure-state  energy  functionals.  The 
problem  of  ensemble  N-representability  versus  pure-state  N-representability  was  also 
discussed  by  different  authors  [64,  67-69].  They  all  came  to  the  conclusion  that  the 
distinction  between  pure-state  N-representability  and  ensemble  N-representability  is 
not  necessary. 

Even  though  the  N-representability  condition  for  1-matrices  is  known,  it  is  not 
quite  clear  what  further  restrictions  must  be  imposed  for  a correct  Euler  equation 


(2.59) 


20 


to  result.  Valone  [67]  gave  an  example  of  an  unacceptable  variation  in  which  the 
partitioning  of  the  primed  and  unprimed  coordinates  was  not  maintained. 

However,  an  orbital  parameterization  of  the  1-matrix  is  always  possible  [25,  63, 
68,  70,  71].  D1(x,  x')  can  be  expanded  in  terms  of  the  natural  orbitals 

OO 

Dl{x,  x')  = ^ np(j>p(x)<f>*{a/).  (2.60) 

p- 1 

This  is  rather  a formal  expression  because  natural  spin  orbitals  are  obtained  a pos- 
teriori. For  this  reason,  some  authors  [63,  68]  expand  the  1-matrix  in  a nondiagonal 
representation  of  orthonormal  spin  orbitals.  However,  if  the  Euler  equations  ob- 
tained with  that  more  general  ansatz  are  cast  into  canonical  form,  the  same  Euler 
equations  are  recovered  as  for  the  diagonal  representation  [63]. 

The  N-representability  conditions  for  1-matrices  in  their  natural  orbital  repre- 
sentation [72]  can  conveniently  be  introduced  with  the  method  of  Lagrange  multi- 
pliers. If  np  = cos2  9P,  the  following  constraints  must  be  applied, 

OO 

(<f>p\<f>g)  = (p\q)  = fipq,  COs2  9P  = N'  0 < < 7T.  (2.61) 

p=  1 

It  is  assumed  that  the  natural  orbitals  form  a complete  set.  The  functional  E{D1} 
is  now  replaced  by 


SiD1}  = E{D1}  - ^2  Kq(q\p)  + A(iV  - ^2nP),  (2.62) 

pq  p=  1 

where  Xpq  and  A are  Lagrangian  multipliers  and  variations  are  with  respect  to  <j>*, 
(j)p  and  0p.  Beginning  with  0*  the  variation  of  the  last  two  terms  of  the  right-hand 
side  of  Eq.  2.62  is  given  by 


S (Er«?  Xrq  f <t>*r(x)<t>q {x)dx} 

Scp*p(x) 

(N-^np)) 

H*p{x) 


^ ] ^pqfiqix) 

Q 


(2.63) 


0. 


(2.64) 
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To  vary  the  first  term  of  the  right-hand  side  of  Eq.  2.62  the  formal  identity 

SE  rr  SE 


= [ [ SE 
J J 5D1(x',x" 


(2.65) 


S(ft*(x)  J J SD1(x',x ")  S(ft*(x) 

can  be  made.  It  should  be  mentioned  at  this  point  that  the  existence  of  8D1(x , x') 

is  not  proven.  However  its  existence  is  a reasonable  assumption  [25]. 

The  variation  of  the  1-matrix  with  respect  to  <ft*  is  given  by 

8D\x',x")  _ 5(E  qMx')nq<l>*q(x")) 

5(p*p(x)  Scft*(x) 

= np(j)p(x,)8(x"  — x). 


(2.66) 

(2.67) 


Substituting  Eq.  2.67  into  Eq.  2.65  leads  to 
8E 


= J J n,Mx')s(x”  - x)dx'dx" 

= / 5D^)^x')dx' 


SDl(x' , x) 

= np  J h(x,  x')(ftp(x')dx', 


with 


h(x,  x')  = 


SE 


(2.68) 

(2.69) 

(2.70) 

(2.71) 


8D1(x',  x) 

The  variation  of  the  functional  2.62  with  respect  to  (j)*p  is  obtained  by  combining 
Eq.  2.70,  Eq.  2.63  and  Eq.  2.64 


SE 


— nph(j)p  \q(t>q  — 0. 


(2.72) 


Similarly,  the  foregoing  procedure  can  be  repeated  for  the  variation  with  respect 
to  (ftp,  leading  to 


se_ 

8 (ftp 


— rip(ftph  ^ ^ pq(ftq  ~ 0. 


(2.73) 


The  two  equations  2.72  and  2.73  can  be  reduced  to  a single  condition 


(rip  - nq){p\h\q)  = 0. 


(2.74) 
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Since  the  energy  is  written  in  terms  of  the  1-matrix,  the  orbitals  in  h occur  in  pairs 
of  (f)*  and  (pp.  Therefore,  the  operator  h is  invariant  under  unitary  transformation 
of  the  orbitals,  and  the  orbitals  can  always  be  chosen  to  be  the  solutions  to  the 
eigenvalue  equation 

— ^p0p>  (2.75) 


where  ep  = (p\h\p). 

Finally  the  functional  2.62  has  to  be  varied  with  respect  to  9P.  The  variation 
of  the  last  two  terms  of  the  right-hand  side  of  Eq.  2.62  is  given  by 


S (Epg  Kq  I <t>*p(x)(t>q(x)dx ) 
56 p 

5 (\(N  — E^lic°s2^r)) 
56p 


0 

A 2 sin  6P  cos  9P 
A sin  29 p. 


(2.76) 

(2.77) 

(2.78) 


Again  the  formal  identity 


Mp  J J 


5E  5D1(x',  x)  . 

-ax  ax 


(2.79) 


5D1(x\x)  59p 

can  be  used,  where  the  variation  of  the  density  matrix  with  respect  to  9V  is  given  by 
5Dx(x',x)  _ 5(£9<MOcos  2090*(x)) 


59n 


59„ 


= -sin  2 9p(pp(x')(f)l(x). 


(2.80) 

(2.81) 


Substitution  of  Eq.  2.81  in  Eq.  2.79  leads  to 


5E_ 
59 p 


sin  29 p II  SD^^Mx'Wp(x)dx,dx 

sin  29pCp, 


(2.82) 

(2.83) 


with 


eP= j 


5E 


5Dx(x',  x) 


4>p{x')dx'  dx. 


(2.84) 
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Combining  Eq.  2.83,  Eq.  2.76,  and  Eq.  2.78  the  variation  of  the  functional 
2.62  with  respect  to  6P  is  obtained 


From  Eq. 
satisfied 


5£ 

— = sin  2 6p(ep  - A)  = 0. 

ot/„ 


(2.85) 


2.85  it  can  be  deduced  that  one  of  the  following  conditions  has  to  be 


(i)  np  = 1 (ii)  ep  = X (iii)  np  = 0, 


(2.86) 


implying  that  all  partially  filled  natural  orbitals  must  belong  to  the  same  degen- 
erate eigenvalue.  There  is  a many-to-one  mapping  between  the  eigenvalues  of  D1 
and  the  eigenvalues  of  h for  partially  filled  orbitals.  These  results  have  been  verified 
by  Donelly  and  Parr  [63]  and  Zumbach  and  Maschke  [70].  If  one  attempts  to  ob- 
tain the  exact  density  matrix  that  is  built  from  partially  occupied  natural  orbitals, 
the  operator  h must  be  quite  different  from  the  one-electron  operators  employed  in 
Hartree-Fock,  and  similar  theories  and  the  equations  must  be  of  radically  different 
mathematical  structure  [73].  The  problem  disappears  in  the  Hartree-Fock  approxi- 
mation for  a closed-shell  system,  because  all  orbitals  are  then  either  fully  occupied 
or  empty.  Even  if  general  one-matrices  are  allowed,  the  solution  of  the  Hartree-Fock 
equations  at  stationarity  is  an  idempotent  density  matrix  [74], 

2.2.4  Hartree-Fock  Theory:  a Special  Case 

Assuming  that  the  correlation  energy  is  zero,  the  expression  for  the  operator  h 
in  Eq.  2.75  can  be  explicitly  found  [71].  The  eigenvalue  equation  was  given  in  Eq. 
2.75 


£p4*p^ 


(2.87) 


with 


5E 

8Dx(x' , x ) 


h(x,  x')  = 


(2.88) 
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In  an  orthonormal  basis  {tpp}  (which  could  for  instance  be  the  basis  of  atomic 
orbitals,  orthogonalized  with  Lowdin’s  procedure  ipp  = SpPXr,  Spr  being  an 
element  of  the  overlap  matrix)  the  matrix  representation  of  the  1-matrix  is  given  by 


PpqVp(x)ip*q(x'),  (2.89) 

v 1 

with  Ppq  being  the  elements  of  a Hermitian  matrix  Ppq  — (p\Dl\q).  Once  the  basis 
has  been  introduced,  the  1-matrix  is  fully  specified  by  the  matrix  P,  and  the  energy 
functional  E{D1}  is  converted  to  a function  of  P.  Then,  the  functional  derivative 
of  E{D1}  with  respect  to  Dl  can  be  written  as  the  partial  derivative  of  the  energy 
with  respect  to  the  elements  of  P 


Hpq<pp{x)ip*q(x'),  (2.90) 

V Q 

with 

Hpq  = ( p\h\q ),  {p\h\q)  = ■ (2.91.) 

If  the  correlation  energy  is  zero,  the  energy  functional  of  the  density  matrix  can 
be  written  as 


E ~ J V2D1(xi,x'1)\Xl=x'idx1  + J v(xi )D1(x1,x1)dx1 

^ J J ^D1(x1,xl)D1(x2,x2)dxidx2 
-J  J ^-D1(x1,x2)D1(x2,xi)dxidx2Sj  . 

In  matrix  representation  the  energy  is  given  by 

E = J2Prs  ((r|-^V2  + n(x)|s)  + ^^Pw(rp||s9)]  . 

rs  \ pq  / 


(2.92) 


(2.93) 
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The  matrix  elements  of  the  matrix  H can  be  found  by  taking  the  derivative  of  the 
energy  with  respect  to  Ptu 

d((r|  - ^V2  + v(x)\s)  + ^pqPpq{rp\\sq)) 

+ 5>'“ 3K ' 

rs 

= {t\  ~ -\-v(x)\u) + ^2Ppq(tp\\uq).  (2.94) 

pq 

The  matrix  H represents  the  nonlocal  Fock  operator,  defined  in  Eq.  2.7. 

2.3  Correlated  Independent  Particle  Model 
The  previous  section  establishes  the  foundation  of  density  matrix  functional 
theory  in  the  sense  that  it  is  possible  to  find  an  expression  for  the  energy  in  terms  of 
the  one-matrix  only,  and  that  the  variational  principle  can  be  applied.  However,  if 
the  one-matrix  is  expressed  in  terms  of  natural  orbitals,  the  resulting  Euler  equations 
reveal  that  all  eigenvalues  of  D 1 that  lie  between  0 and  1,  must  be  mapped  into  a 
single  eigenvalue  of  the  one-electron  operator  h(x,x').  To  avoid  the  degeneracy 
problem,  the  density  matrix  within  the  CIP  model  is  restricted  to  be  idempotent. 
The  operator  h(x,x')  includes  a correlation  potential  vc{Dx)  = so  that 

the  orbitals  and  the  eigenvalues  of  the  CIP  model  include  correlation.  The  density 
matrix  cannot  be  expected  to  be  exact.  However,  the  question  of  how  good  the 
density  matrix  obtained  from  a single  determinant  can  be  will  be  addressed  in  a 
later  chapter. 

If  the  correlation  energy  Ec  is  zero,  the  energy  functional  E{D1}  is  given  by 
the  Hartree-Fock  energy  functional.  If  the  correlation  energy  is  nonzero,  and  the 
density  matrix  restricted  to  be  idempotent,  the  operator  h(x,  x')  can  be  expressed 
as  the  sum  of  the  Fock  operator  f(x,x')  and  an  operator  vc(x,  x'),  the  correlation 
potential 


h(x,  x')  = f(x,  x')  + vc(x,  x'). 


(2.95) 
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If  the  correlation  potential  is  derived  from  a partitioned  equation-of-motion 
approach,  as  described  in  the  next  chapter,  the  correlation  potential  is  energy- 
dependent.  The  resultant  orbitals  are  nonorthogonal,  and  the  density  matrix  is 
therefore  not  idempotent.  Strictly  speaking,  the  model  is  not  an  independent  particle 
model,  if  an  energy-dependent  correlation  potential  is  used,  and  the  density  matrix 
could,  in  principle,  be  obtained  exactly.  The  disadvantage  of  an  energy-dependent 
correlation  potential  is  that  a different  eigenvalue  problem  has  to  be  solved  for  each 
orbital  separately  with  corresponding  eigenvalue.  Our  efforts  therefore  aim  towards 
an  energy-independent  correlation  potential,  which  can  be  constructed  using  the 
Fock  space  coupled-cluster  method. 

Our  goal  is  to  express  the  correlation  potential  vc  in  a basis  (i.e.,  a matrix 
representation),  symbolized  by 

Vc  = H - F.  (2.96) 

Since  the  eigenvalues  of  F are  Koopmans’  values  for  ionization  potentials  (IP’s)  and 
electron  affinities  (EA’s),  a reasonable  assumption  would  be  that  the  eigenvalues  of 
H are  the  exact  IP’s  and  EA’s.  The  condition  of  exactness  of  the  IP’s  and  EA’s  as 
the  eigenvalues  of  H could  then  be  used  to  derive  the  correlation  potentials  for  the 
CIP  model. 

The  extended  Koopmans’  theorem  employs  one-particle  operators  whose  eigen- 
values are  IP’s  and  EA’s  while  using  the  exact  wave-function.  Since  the  extended 
Koopmans’  theorem  (EKT)  inspired  and  supported  our  choice  of  the  condition  we 
want  to  impose  on  the  correlation  potential  of  the  CIP  model,  the  EKT  is  described 
next. 

2.4  Extended  Koopmans’  Theorem 

Within  the  Hartree-Fock  theory,  Koopmans’  theorem  Eq.  2.10  and  Eq.  2.11 
provides  a simple  one-electron  model  for  ionization  and  electron  attachment.  A 
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similar  model  using  the  exact  wave-function  or  a correlated  wave-function  for  the 
unionized/unattached  system  is  given  by  the  extended  Koopmans’  theorem,  which 
was  independently  developed  by  Day,  Smith  and  Garrod  [52,  59,  75,  76]  for  single 
ionization  and  attachment  processes  and  a generalization  of  it  for  ionization,  which 
includes  double  and  higher  ionization  by  Morrell,  Parr  and  Levy  [53]. 

The  EKT  method  folds  the  information  of  the  2-matrix  into  two  distinct  one- 
particle  potentials.  One  has  eigenvalues  that  represent  ionization  energies  and  eigen- 
functions, which  represent  the  orbitals  from  which  the  electron  is  removed.  The  other 
has  eigenvalues  that  represent  electron  affinities  and  eigenfunctions,  which  represent 
the  orbitals  for  the  addition  of  an  electron. 

The  basic  assumption  of  the  EKT  method  is  that  ionization  and  electron  capture 
can  be  described  by  the  creation  or  annihilation  of  a single  spin  orbital.  Ionization 
is  viewed  as  elimination  of  a variationally  determined  orbital  from  a fixed  correlated 
reference  wave-function  of  the  N-electron  system,  electron  attachment  as  addition 
of  a variationally  determined  orbital  to  a fixed  correlated  reference  wave-function. 
The  model  accounts  for  correlation  and  relaxation  effects. 

The  wave-function  of  the  N + 1 system,  where  an  electron  in  orbital  ipr  was 
added  to  the  wave-function  of  the  N electron  system  can  be  expressed  as 

^+1  = (N  + 1)~^AM^N(X2  ■ ■ ■ xN+i)  = 0^N,  (2.97) 

where  r is  the  particular  state  of  the  negative  ion.  The  wave-function  of  the  N — 1 
system,  where  an  electron  in  orbital  ips  was  deleted  from  the  wave-function  of  the 
N electron  system,  is  given  by 

= ^ J ra(x1)VN(x1,x2...xN)dx1=OlVN , 


(2.98) 
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where  s is  the  particular  state  of  the  positive  ion  and  A is  the  antisymmetrizer.  The 
operators  O*  and  the  optimized  one-electron  function  are  defined  as 


o\ 

= Ylpjctr 

P 

ot 

= 

V 

(2.99) 

4>r 

P 

i>s 

V 

(2.100) 

where  p*  is  the  creation  operator  and  p is  the  annihilation  operator  for  the  basis 
spin  orbital  <pp)  and  the  cps  and  the  c+r  are  the  coefficients,  which  are  to  be  deter- 
mined variationally.  The  creation  and  annihilation  operators  satisfy  the  following 
anticommutation  relations 


0 = sr  -\-rs 

0 = sM  + 

Srs  = sr t + r^s. 


(2.101) 


If  TiV  is  normalized,  and  u is  r or  s,  then  the  energy  difference  between  the 
ionic  and  the  parent  species  is  defined  by 


(OlVN\H\OlVN) 

(olvN\olvN) 


(Vn\H\Vn). 


(2.102) 


is  an  orbital-energy-like  quantity,  because  if  T;V  is  frozen  and  t/t  is  varied,  the 
right-hand  side  of  Eq.  2.102  depends  only  on  the  parameters  in  the  orbital  Eq. 
2.102  takes  a simple  form,  if  TiV  is  an  eigenfunction  of  the  Hamiltonian  H,  because 
then  the  second  term  on  the  right-hand  side  of  Eq.  2.102  can  be  written  as 


n1h^N)  = 

{tyN\OuOl\VN){VN\VN)E 

(olvN\ol\yN) 

(*N\OuOtH\*N) 

{ol^N\ol\^N)  ’ 


(2.103) 
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and  Eq.  2.102  becomes 


(2.104) 


In  case  of  an  electron  removal,  O £ = 0\,  given  by  Eq.  2.99.  Insertion  into  Eq. 
2.104  leads  to 


(2.io5) 

pq  pq 

which  can  be  written  in  matrix  form 

(c-^V-  - e"S-)c-  = 0,  (2.106) 

with  the  matrices  V~  and  S“  defined  as 

V'  = -<*'V[tf,«]l*JV>  (2.107) 

5-  = (2.108) 

v-  is  the  effective  one-body  potential  in  which  the  electron  to  be  ionized  moves. 
The  variationally  stable  orbitals  c“  are  the  eigenvectors  of  V-  with  respect  to  the 
metric  S~.  The  eigenvectors  and  eigenenergies  are  solutions  to  the  equation 

Vc_  = S~c-e-.  (2.109) 

The  matrix  S'  can  be  identified  with  the  first-order  reduced-density  matrix 
written  in  second  quantization.  With  the  metric  S~  the  domain  of  the  operator  V” 
is  limited  to  the  region  of  the  one-particle  space  spanned  by  the  eigenfunctions  of 
the  1-matrix,  having  nonzero  eigenvalues. 

If  the  reference  function  is  not  an  exact  eigenfunction  of  the  Hamiltonian,  V~ 
will  not,  in  general,  be  Hermitian.  The  right-hand  side  of  Eq.  2.107  should  then  be 
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replaced  by  its  Hermitian  component 

V~  = + [p*.  (2-110) 

The  metric  remains  unaltered. 

To  derive  an  explicit  expression  for  V-,  the  Hamiltonian  in  second  quantization, 
given  by 


H 


^ hrsr^S  + 

rs 

hrsr^s  + 


1 

2 

1 

4 


y ^(rs\ut)r^  sHu 

rstu 


(2.111) 


rs  rstu 

is  inserted  in  Eq.  2.107.  Using  the  anticommutation  relations  2.101,  the  one-electron 
part  of  V~  can  be  written  as 


,N\ 


KT  = 

rs 

= - ^ hrs('l!N\p*rhq\'f'N)  + ^ hrs(TAf|ptqrts|4'iV) 

rs  rs 

= Y,  Ka(VN\p*s\'J'N)Sq 

rs  rs 

- yy  M^b^bsi^) = yy  hqsDivs. 

rs  s 

Similar  steps  can  be  done  for  the  two-electron  part 


Jqr 


(2.112) 


K)  = -^y^(^i«o<^jvbt[?'tst^9]i^yv) 

rstu 

= — - ^^(rs\ut)('f>N\p^r^ shuq\tyN)  + - y~]  (rs|ut)  (TjV[pbrb^u|^>iV) 
z z 

rstu  rstu 

= — - yy(rs|ut)(T7V|ptrbt^ul^Ar)  H — y^(rs|uf)(TjV|pb^u|TAr)(5gr 
2 2 

rstu  rstu 

— - yy(rs|ut)(^,ivipV^u|^/v)59S  + - yy(rs|ut)(^yv|pvtsb^|4,iv) 

2 2 
rstu  rstu 

= = 2y^<r9lb0^tu> 

rtu  rtu 


(2.113) 
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where  D2rtu  is  the  symmetric  part  of  the  second-order  reduced-density  matrix 

- Dlut.  (2.114) 

Adding  Eq.  2.112  and  Eq.  2.113  leads  to 

Vip  = KDls  + 2 J2(rQ\\ut)Dlrtu-  (2-115) 

s rtu 

For  exact  2-matrices  D2rtu  = —Dprut  holds,  for  approximated  wave- functions  that 
might  be  an  approximation.  The  matrix  V is  also  known  as  the  generalized  Fock 
matrix  [66]. 

The  coordinate  representation  of  the  nonlocal  potential  Eq.  2.115  is  given  by 

V~(x,  x')  = h(x)D1(x,  x')  + 2 J u(x,x")D2(x"x,x"x')dx".  (2.116) 

Eigenfunctions  and  eigenenergies  to  that  operator  are  solutions  to 

J V~(x,  x')^s(x')dx'  = es  J Dl(x,  x')ij}a{x').  (2.117) 

The  one-particle  functions  iJjs  are  some  linear  combination  of  natural  orbitals  having 
nonzero  occupancy  in  the  reference  function  T iV . 

Similarly,  for  the  attachment  of  an  electron  the  one-electron  potential  is  defined 

by 

K = (2.118) 

= h„  + ^(P<ll?“)Oi,-V',;.  (2.119) 

tu 

The  potential  is  defined  with  respect  to  the  metric 

S*  = (fl'Wi#'') 

= dqp  ~ Dqp 


(2.120) 

(2.121) 
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The  corresponding  nonlocal  operator  to  Eq.  2.119  in  coordinate  representation  is 
given  by 

V+(x,x')  = 8{x  - x')h(x)  + S(x  - x')  J u(x,  x")Dl{x"  ,x")dx"  - u(x,  x')Dl{x,x') 
— h(x)D1(x,  x')  — 2 J u(x,x")D2(x"x,x"x')dx".  (2.122) 

The  metric  in  coordinate  representation  can  be  written  as 

5+(x,  x')  = 8(x  — x')  — B1(x,  x').  (2.123) 

The  potential  is  dehned  with  respect  to  the  orbital  basis  spanned  by  natural  orbitals, 
having  nonunity  occupation  numbers.  If  the  2-matrix  is  not  derived  from  an  exact 
eigenfunction  of  the  Hamiltonian,  the  Hermitian  component  is  used  as  before. 

Defining  a new  potential  by  combining  the  first  three  terms  of  the  right-hand 
side  of  Eq.  2.122 

V (x,  x1)  = 8(x  - x’)h(x)  + 5(x  - x')  J u(x,  x")Dl(x",  x")dx"  - u(x,  x,)D1(x,  x'), 

(2.124) 

the  one-particle  potential  V+(x,  x')  can  be  written  as 

V+{x , x')  = V (x,  x')  - V~(x , x1).  (2.125) 

V corresponds  directly  to  the  Hartree-Fock  potential  with  the  Hartree-Fock  density 
matrix  replaced  by  a correlated  1-matrix. 

The  variationally  stable  orbitals  to  which  an  electron  is  added,  and  the 
eigenenergies  are  obtained  as  solutions  of 

J V+(x,x')ipr(x')dx'  = er  J S+(x,x')i>r(x '),  (2.126) 

which  in  matrix  representation  is  given  by 


V+c+  = S+c+e+. 


(2.127) 
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Orthogonalization  of  the  eigenvectors  by  substituting  c1*1  = S*  ^(c')*  into  Eq. 
2.109  and  Eq.  2.127,  and  symbolizing  ionization  or  attachment  by  ±,  leads  to 

F IktM*  = (c')±e±  (2.128) 

with 

F±KT  = S±-iv±S±“*.  (2.129) 

fekt  the  nonlocal  extended  Koopmans’  theorem  Fock  operator  in  matrix  form. 
Fekt  is  nonhermitian  for  approximate  wave-functions  that  include  correlation.  How- 
ever, its  eigenvalues  are  real  and  equal  to  the  variationally  optimized  energy  differ- 
ences between  the  ionic  and  parent  species:  the  IP’s  and  the  negative  of  the  EA’s. 

The  question  whether  or  not  the  EKT  method  is  exact  has  been  discussed  by 
several  authors.  Morrell  et  al.  [53]  suggested  that  the  lowest  ionization  potential 
is  correct.  Katriel  and  Davidson  [54]  showed  that  if  the  N-electron  wave-function 
is  exact,  at  least  the  smallest  ionization  energy  and  the  corresponding  orbital  are 
exact,  meaning  that  the  partner  to  the  eigenfunction  of  the  one-particle  operator 
is  an  exact  eigenfunction  of  the  (N  — l)-particle  Hamiltonian.  Based  on  a second- 
order  perturbation  analysis  Pickup  and  Snijders  [77]  concluded  that  except  for  a 
two-electron  system  none  of  the  EKT  eigenvalues  is  exact.  Numerical  support  for 
the  exactness  of  the  lowest  ionization  potential  was  given  by  Morrison  and  Sundholm 
and  Olsen  [56-58] . Olsen  and  Sundholm  [55]  also  argued  that  for  the  lowest  IP  and  a 
complete  basis  set,  the  second  and  higher-order  correction,  obtained  from  the  EKT 
method,  are  correct.  Recently  Pernal  and  Cioslowski  [78]  published  a theorem  that 
provides  a sufficient  condition  for  the  validity  of  the  extended  Koopmans’  theorem. 

The  EKT  defines  one-particle  operators,  whose  eigenvalues  are  IP’s  and  the 
negative  of  EA’s,  of  which  at  least  the  highest  IP  is  exact,  if  the  exact  wave-function 
is  used.  That  supports  our  choice  of  the  condition  on  the  correlation  potential  to 
reproduce  the  exact  IP’s  and  EA’s  as  eigenvalues  of  h{xi,x\). 
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If  the  wave-function  is  a single  Slater  determinant,  the  2-matrix  can  be  written 
as 


D2(xxx2,x\x'2)  = [D1(xi,x'1)D1(x2,x'2)  - D1(xi,x,2)D1(x2,x'1)\  . (2.130) 

The  one-particle  potentials  reduce  to 

K-(X.X')  = Jv(X,XW,X')dX"  (2.131) 

V+M  = (2.132) 

These  are  representations  of  the  Hartree-Fock  potential,  projected  onto  the  disjoint 
spaces  spanned  by  D1  and  1 — D1.  Orthonormalization  of  the  eigenvectors  gives 
the  Fock  operator  in  its  common  form.  The  eigenenergies  of  the  two  potentials  are 
Koopmans’  values  for  IP’s  and  EA’s. 

2.5  Alternative  Models  Based  on  the  Hohenberg-Kohn-Gilbert 

Theorem 

Within  the  last  six  years,  density  matrix  functional  theory  based  on  the  Hohenberg- 
Kohn-Gilbert  theorem  has  received  new  interest.  Most  of  the  proposed  density  ma- 
trix functionals  are  written  in  terms  of  natural  orbitals,  and  are  usually  referred  to 
as  natural  orbital  functionals. 

The  total  energy  can  be  expressed  in  terms  of  the  natural  orbitals  (j)p,  and  the 
diagonal  elements  of  the  second-order  reduced-density  matrix  er(xi,  x2) 


a(xi,x2)  = D2(x1x2,xix2). 


(2.133) 


The  energy  is  given  by 

E — J 0p(^i)V20p(xi)dxi  + J v(xi)D1(xi,xi)dxi 


+ 


11 


cr(xux2) 


r 12 


dx\dx2. 


(2.134) 
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In  order  to  formulate  a functional  that  depends  on  natural  orbitals,  the  two-particle 
density  a has  to  be  approximated,  since  the  exact  dependence  of  the  2-matrix  on 
the  1-matrix  is  not  known.  The  different  natural  orbital  functionals  differ  in  the 
approximation  scheme  used  for  the  two-particle  density. 

Csanyi  and  Arias  [26]  introduced  the  terminology  of  corrected  Hartree  and 
corrected  Hartree-Fock  functionals,  in  which  most  of  the  natural  orbital  functionals 
can  be  divided.  They  considered  tensor  product  approximations  to  the  2-matrix  to 
provide  estimates  of  the  2-matrix  in  terms  of  the  1-matrix.  The  2-matrix,  which 
depends  on  four  variables,  is  expanded  in  terms  of  one-body  operators  gp  and  hp, 
which  depend  on  two  variables, 

D2  = ^2gp®hp.  (2.135) 

p 

The  tensor  product  denotes  one  of  the  three  possible  choices  for  separating  four 
variables 


7 : 

g ® h = g(xi,x[)h(x 2,x2) 

(2.136) 

77  : 

g®h  = g(xi,x2)h(x'1,x'2 ) 

(2.137) 

III  : 

g <S>  h = g(x1,x'2)h(x2,x[). 

(2.138) 

If  the  2-matrix  is  approximated  by  a single  tensor  product  of  type  7,  the  familiar 
Hartree  approximation  is  recovered,  and  a corrected  Hartree  scheme  can  be  defined 
by 

dch  = 5 (D'  ®7  D')  + Dle-  (2.139) 

The  factor  | is  needed  to  maintain  normalization.  D2C  represents  exchange  and 
correlation  effects,  and  can  be  expanded  according  to  Eq.  2.135. 
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Alternatively  antisymmetry  might  be  ensured  explicitly,  which  leads  to  the 
corrected  Hartree-Fock  scheme 

Dchf  ~ 2 (■^1  — D 1 (&ui  D1)  + D2.  (2.140) 

The  unknown  D2,  representing  the  correlation  effects,  can  again  be  expanded  as  in 
Eq.  2.135. 

Terms  of  type  II  violate  the  condition  that  D2  — > 0 as  the  insertion  x\x'2 
and  removal  X1X2  locations  are  placed  further  apart.  Such  terms  are  therefore  not 
included  in  an  expansion  of  D 2 . 

If  one  wants  to  expand  D 2 beyond  the  Hartree-Fock  approximation,  antisym- 
metry requires  representing  D2  as  a pair  of  terms,  so  that  in  its  simplest  form 

1)2  = \ ( D*  ®/  Dl  ~ Dl  ®///  Dl  + 9 ®i  9 ~ 9 ®///  g)  ■ (2.141) 

The  sum  rule  requires  g — y/ ±DX(  1 — D1).  However,  the  structure  of  the  energy 
functional  Eq.  2.134  is  such,  that  for  this  choice  the  resulting  energy  functional  of 
D 1 is  equal  to  that,  generated  by  representing  D2  as  a single  type  III  product. 

The  two  simplest  representations  of  the  2-matrix  to  explore  are  the  corrected 
Hartree  approximation 

Dch  = 2 ~ V D1  <8>ni  Vm)  (2.142) 

and  the  corrected  Hartree-Fock  approximation 

D2chf  = l (D1  ®7  D1  - D1  ®ni  D1  - y/D'{l-D')  ®IU  - D1))  . 

(2.143) 

Expressing  the  density  matrix  in  terms  of  natural  orbitals,  leads  to  the  different 
natural  orbital  functionals. 

However,  there  are  different  choices  for  the  expansion  of  the  2-matrix  possible, 
if  more  terms  are  included  in  the  expansion  Eq.  2.135.  Holas  [79]  also  suggested 
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expanding  the  2-matrix  as  a linear  combination  of  products  of  type  I and  III,  which 
leads  to  the  corrected  Hartree  and  corrected  Hartree-Fock  approaches  as  special  cases 
of  such  an  expansion. 

Godecker  and  Umrigar  [28]  employed  a corrected  Hartree  type  expression  for  the 
2-matrix.  For  a spin-independent  Hamiltonian  each  natural  orbital  can  be  chosen 
to  be  either  spin  up  or  spin  down,  and  be  labeled  by  an  orbital  index  p and  a spin 
index  sp.  The  approximate  a in  Eq.  2.134  takes  the  following  form 


a[{n }>  W]  = Y ~ Y H q SsPsMri)Mri)Mr2)Mr2)- 


PQ 


pq 


(2.144) 

The  np  are  the  occupation  numbers.  The  primes  indicate  that  the  p — q terms  are 
omitted.  Omitting  the  terms  p = q in  the  sums  in  Eq.  2.144,  introduces  a correction 
for  electron  self-interaction,  but  it  also  causes  the  sum  rule  not  to  be  fulfilled. 

The  ground  state  can  be  found  by  substituting  the  definition  of  a , Eq.  2.144, 
into  Eq.  2.134  and  minimizing  the  energy  with  respect  to  the  natural  orbitals  and 
the  occupation  numbers  under  the  constraint  that  the  natural  orbitals  be  orthogonal. 
The  functional  derivatives  are 


dE 

<%(o) 


dE 

Orir, 


u 7 r 

= V20P(o) + npu(ri)0p(r!)  + ^npn90p(ri)  / 

Q 

~ ^ ^ \/npnq^svsq4>q{r\)  J 


02 


0p(o)0g(r2) 


drq 


(2.145) 


= \ J + J w(ri)0j(o)dri 


^(o)^(o) 


-t Urj [ J 

q J J 


dr i dr 2 


0p(o)0q(o)0p(ri)09  (r ! ) 


02 


dodr2.  (2.146) 


Note  that  the  resulting  functional  coincides  with  the  Hartree-Fock  functional,  if  the 
occupation  numbers  are  constrained  to  be  0 or  1. 
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Buijse  and  Barends  [27,  80]  independently  derived  an  expression  for  the  two- 
particle  density  that  is  equivalent  to  Eq.  2.144,  but  where  the  sums  run  over  all 
indices  i and  j,  and  is,  therefore,  not  electron  self-interaction  corrected.  Eq.  2.144 
was  derived  by  approximating  the  full  exchange  and  coulomb  correlation  hole  in  a 
way  that  is  analogous  to  the  Fermi  hole  in  the  Hartree-Fock  model. 

Despite  their  simplicity,  the  corrected  Hartree  (CH)  and  corrected  Hartree- 
Fock  (CHF)  functionals  predict  reasonably  accurate  correlation  energies  for  atoms 
[28,  81].  On  the  other  hand,  the  reconstructed  2-matrices  fulfill  only  a few  necessary 
conditions  for  N-represent ability  [32,  33,  79],  and  fail  to  describe  the  homogeneous 
electron  gas  properly[29-31].  It  was  reported  that  the  CH,  CHF,  and  hybrids  [32] 
potential  energy  curves  for  diatomics  have  almost  no  minimum,  because  of  severe 
overcorrelation  at  large  R [32,  82],  whereas  the  Godecker  and  Umrigar  functional 
gave  plausible  molecular  dissociation  curves  [32]. 

An  alternative  idea  to  construct  the  correlation  energy  functional  was  presented 
by  Yasuda  [34].  He  used  the  density  equation  [83,  84],  whose  basic  variable  is  the 
reduced-density  matrix  instead  of  the  wave-function.  The  density  equation  is  equiva- 
lent to  the  Schrodinger  equation.  Approximation  of  the  correlation  energy  functional 
is  introduced  as  the  decoupling  of  the  higher-order  reduced-density  matrices.  Ya- 
suda proposed  a two-step  procedure.  In  the  first  step  the  first-order  density  equation 
is  used  to  reconstruct  the  one-body  operator  of  the  Hamiltonian.  Once  the  Hamil- 
tonian is  obtained,  the  connected  piece  of  the  second-order  reduced-density  matrix 
A2  is  calculated  by  the  second-order  density  equation,  although  it  could  be  obtained 
by  any  other  ab  initio  method.  The  two  steps  are  repeated  until  convergence  of  A2 
is  reached.  The  correlation  energy  can  then  be  calculated  from  the  converged  A2. 
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The  first-order  reduced-density  equation  is  given  by 

0 = {h(xi)  - E}  Dl(x[,xi) + 2 J {h(x2) +u(xi,x2)}  D2(x'lx2,xix2)dx2 
+ 3 J u(x2,X3)D3(x'1X2X3,xiX2X3)dx2dx3.  (2.147) 

The  reduced-density  matrices  of  the  exact  eigenstate  of  H satisfy  this  equation.  The 
second-  and  third-order  reduced-density  matrices  can  symbolically  be  expressed  as 

D2  = D1  A D1  + A2  (2.148) 

D3  = D1  A Dl  A D1  + 3D1  A A2  + A3.  (2.149) 

The  wedge  product  A generates  the  normalized  antisymmetrized  product.  A3  is  the 
connected  piece  of  the  third-order  reduced-density  matrix.  In  Yasuda’s  approach  the 
third-order  reduced-density  matrix  is  approximated  in  terms  of  D 1 and  A2  alone, 
when  substituted  in  Eq.  2.147. 

If  the  natural  spin  orbitals  are  used  as  the  one-particle  basis,  the  diagonal 
elements  of  the  generalized  Fock  operator  are  given  by 

with  rii  being  the  natural  occupation  numbers.  The  rescaling  factor  a is  introduced 
to  compensate  for  the  effect  of  the  higher-order  terms. 

Eq.  2.147  and  Eq.  2.150  determine  the  one-body  operator  or  the  generalized 
Fock  operator  in  terms  of  the  1-matrix  and  a trial  2-matrix.  With  the  reconstructed 
Hamiltonian,  the  2-matrix  can  be  determined  and  the  calculation  of  the  one-body 
operator  is  repeated  until  convergence.  The  correlation  energy  is  given  by 

Ec  = J u(xi,  x2)A2(xiX2,  xix2)dxidx2 

\ > Tip  (1  Tip ) €p 

2a(2np  — 1) 


(2.151) 
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If  the  second-order  density  equation  is  used  to  determine  the  2-matrix,  then,  in 
principle,  the  exact  correlation  functional  could  be  obtained,  if  one  were  to  use 
the  exact  reconstruction  formulas  of  the  third-  and  fourth-order  reduced-density 
matrices. 

The  Yasuda  functional  reproduces  the  correlation  energy  of  selected  atoms  and 
molecules  accurately  [34],  satisfies  the  homogeneous  scaling  relation  [85],  and  the 
particle- hole  symmetry  [34].  It  also  performs  well  for  the  homogeneous  electron 
gas  [86].  However,  if  the  approximate  reconstruction  formulas  for  the  higher-order 
density  matrices  are  used,  the  N-representability  condition  on  the  2-matrix  is  not 
fulfilled  [34],  and  the  energy  functional  becomes  unbounded  from  below  [87].  Yasuda 
also  presented  a local  approximation  to  the  correlation  energy  functional  [36] , which 
is  very  similar  in  appearance  to  the  local  density  approximation  [3,  88]  in  DFT. 

Mazziotti  also  suggested  an  iterative  scheme  in  which  the  first-order  density 
equation  is  employed  to  construct  an  energy  functional  of  the  1-matrix  [35].  He 
proposed  the  use  of  an  antisymmetrized  geminal  power  ansatz  [61]  to  reconstruct 
the  reduced-density  matrices  of  different  order. 

Alternatively,  Cioslowski  and  Lopez-Boada  [38,  39]  presented  a scheme  that 
applies  the  hypervirial  theorem,  yielding  a functional  of  the  1-matrix.  The  functional 
is  parameterized  by  a single  screening  function,  retrieved  from  the  correlation  energy 
of  the  homogeneous  electron  gas. 


CHAPTER  3 

CORRELATION  POTENTIALS  YIELDING  EXACT  IONIZATION 
POTENTIALS  AND  ELECTRON  AFFINITIES 

In  the  last  chapter  the  correlated  independent  particle  (CIP)  model,  which 
can  be  viewed  as  a generalization  to  the  Hartree-Fock  model,  or  an  approximation 
to  density  matrix  functional  theory,  was  discussed.  Within  the  CIP  model  the 
exchange  energy  and  the  kinetic  energy  can  be  expressed  exactly  in  terms  of  the 
density  matrix.  The  remaining  unknown  is  the  functional  of  the  correlation  energy 
in  terms  of  the  density  matrix.  A model  was  proposed  in  which  the  energy  functional 
is  decomposed  into  the  Hartree-Fock  functional,  plus  a correlation  potential.  From 
Koopmans’  theorem  and  the  extended  Koopmans’  theorem  it  was  found  reasonable 
to  impose  the  condition  of  exactness  of  the  ionization  potentials  (IP)  and  electron 
affinities  (EA)  as  the  negative  of  the  orbital  energies,  in  order  to  obtain  a correlation 
potential. 

In  this  chapter,  I will  discuss  how  to  obtain  a correlation  potential,  yielding  the 
exact  IP’s  and  EA’s.  First  the  equation-of-motion  coupled-cluster  (EOM-CC)  ap- 
proach is  reformulated  to  extract  an  energy-dependent  correlation  potential.  Then 
the  equivalence  to  the  propagator  approach  will  be  shown.  Finally  the  Fock  space 
coupled-cluster  (FSCC)  method  is  exploited  to  obtain  an  energy-independent  cor- 
relation potential  that  is  universal  for  all  orbitals. 

3.1  Correlation  Potential  Derived  from  the  Equation-of-Motion 

Coupled-Cluster  Method 

The  EOM-CC  method  gives  the  exact  IP’s  and  EA’s,  if  no  truncation  is  em- 
ployed. In  order  to  find  a correlation  potential  for  the  CIP  model,  so  that  the 
negative  of  the  independent  particle  eigenvalues  are  the  exact  IP’s  and  EA’s,  the 
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EOM-CC  method  is  reformulated  in  this  section.  An  effective  Hamiltonian  is  ob- 
tained, from  which  a correlation  potential  can  be  extracted.  First  the  coupled-cluster 
(CC)  method  is  introduced,  followed  by  the  EOM-CC  method.  The  partitioning 
technique,  which  yields  the  effective  Hamiltonian,  from  which  an  energy-dependent 
correlation  potential  is  obtained,  is  discussed,  and  finally  the  connection  to  the 
propagator  method  is  shown. 

3.1.1  Coupled-Cluster  Method 

Before  outlining  the  CC  method,  the  concept  of  normal  ordering,  which  is  con- 
veniently used  later,  is  introduced.  A product  of  creation  and  annihilation  operators 
is  said  to  be  in  normal  order  relative  to  the  Fermi  vacuum  if  all  particle  creation 
operators  oA  b\  . . . and  all  hole  annihilation  operators  i,j , . . . are  to  the  left  of  all 

particle  annihilation  operators  a,b,...  and  of  all  hole  creation  operators  i\j\ 

The  Fermi  vacuum  expectation  value  of  a normal- ordered  product  of  such  operators 
vanishes. 

The  Hamiltonian  in  second  quantization  was  given  in  Eq.  2.111.  The  Hamil- 
tonian is  rewritten  in  normal  order  with  respect  to  the  Fermi  vacuum  by  applying 
Wick’s  theorem  (for  a detailed  proof  [89]),  which  states,  that  the  normal-ordered 
result  of  the  product  of  two  normal-ordered  operator  strings  {A}  and  {B},  equals 
the  sum  of  the  normal  product  {AB}  and  all  the  possible  contractions  between 
operators  from  A and  operators  from  B.  The  contraction  of  two  fermion  construc- 
tion operators  is  defined  as  their  anticommutator,  if  they  are  not  in  normal  order 
and  zero  otherwise.  The  fermion  construction  operators  fulfill  the  anticommutation 
relations,  Eq.  2.101. 

A contraction  is  denoted  by  a bullet  following  the  operator,  if  necessary  a double 
bullet  is  used  to  distinguish  between  two  different  contractions.  The  one-electron 
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part  of  the  Hamiltonian  then  reads 


^ ] hrsr^s  — ^ ] hrs^v  s}  4"  ^ h,rsv^  s — ^ ']  4-  ^ ] ha. 

rs  rs  rs  rs  i 

The  two-electron  part  can  be  rewritten  as 

-^^(rs\ut)r^sHu  = — ^^(rs\ut){r^  shu}  H — ^^(rs\ut)  (r^sh'u 


(3-1) 


rstu 


rstu  rstu 

4-  r^shu*  + rtst*t*w  + rV’iu’  + r^s^fu" 


+ 


4- 


rW'u*) 

- ^^(rs|ut){rtsttw}  4-  ((ir|is)  — (ir|si))  {r^s} 

rstu  rsi 

sXlMv'lv)  - <y 


(3.2) 


V 


The  energy  of  the  Hartree-Fock  reference  state  is  defined  as 

Bo  = (4>o|tf  l®„>  = Z fc»  + \ E «y  I w>  - <v li«»  . (3-3) 

i ij 

and  the  matrix  elements  of  the  Fock  matrix  corresponding  to  To  are  given  by 

frs  = firs  + ((ir|zs)  - (ir|si)) . (3.4) 

i 

Combining  Eq.  3.1,  Eq.  3.2,  Eq.  3.3,  and  Eq.  3.4,  the  Hamiltonian  is  rewritten  as 


H = E0  + ^2frsWs}  + ^^2(rs\ut){r]shu} 

rs  rstu 

= E0  + ^frsi^s}  + ^ ^{rs||nt){rtsttu},  (3.5) 

rs  rstu 

The  normal-ordered  Hamiltonian  H ^ can  now  be  defined  as 

Hn  — H - Eq  — ^2frs{r^s}  + ^^(rsHut^rVtii}.  (3.6) 

rs  rstu 

Having  the  necessary  conventions  established,  the  coupled-cluster  ansatz  can 
be  discussed.  The  coupled-cluster  wave-function  [90,  91]  is  parameterized  by  an 
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exponential  ansatz 

^ = eT$0,  (3.7) 

where  T is  an  excitation  operator.  The  exponential  ansatz  ensures  extensivity,  and 
exhibits  better  convergence  to  the  full  configuration  interaction  (Cl)  limit  than  a 
linear  expansion. 

The  operator  T consists  of  one-body  (Ti),  two-body  (T2),  etc.  cluster  operators 

T = Ti+T2  + T3 (3.8) 

Each  cluster  operator  is  a sum  of  creation  and  annihilation  operators,  which  fulfill 
Eq.  2.101,  multiplied  by  an  appropriate  coefficient,  cailed  the  amplitudes 

Tl  = £ <?{«'!} 

ia 

ijab 

T3  = 4 £ tgE  {oWjc't}  ....  (3.9) 

ijkabc 

The  braces  {}  indicate  normal  ordering.  In  general  the  cluster  operator  is  given  by 

~ (rn!p  ^ ^...  • • •}>  (3.10) 

ij  - 
ab... 

where  m < N,  the  number  of  electrons.  The  factor  t-ttt  accounts  for  the  fact,  that 
permutations  among  the  hole  indices  or  the  particle  indices  do  not  create  distinct 
contributions.  For  example 

j.ab  j.ab  j.ba  j.ba 

Lij  ~ lji  ~ ~ 

aUtfj  = —a)jbH  — —bHo)j  = b^  jo)  i (3-11) 


contribute  4 equal  terms. 
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If  the  operator  T is  not  truncated,  the  expansion  Eq.  3.7  is  the  exact  wave- 
function.  However,  to  make  the  calculation  feasible,  the  operator  T usually  has  to 
be  truncated.  A common  approximation  is  the  CCSD  method,  in  which  only  Ti  and 
T2  are  kept  in  Eq.  3.7.  If  the  exponential  is  expanded  in  a Taylor  series 

eT  = 1+T  + ^T2  + — T3  + ...,  (3.12) 

Z u! 

the  CCSD  wave-function  is  given  by 

^ = $o  + Ti$o  + o + -Tj2$0  + 7\T2$o  + 2^*0 

+ + \rlm 0 + \t,T^ o + ^T23$o  + . . . . (3.13) 

Terms  in  T,  which  contain  only  a single  cluster  operator,  are  called  connected- 
cluster  contributions.  Terms,  which  contain  products  of  cluster  operators,  are  called 
disconnected  contributions,  which  are  the  reason  for  the  extensivity  of  the  coupled- 
cluster  method. 

The  coupled-cluster  equations  can  be  derived  with  the  help  of  Wick’s  theorem. 
The  general  idea  holds  for  any  truncation  level,  however,  the  following  material 
will  concentrate  on  the  CCSD  approximation.  Inserting  the  coupled-cluster  wave- 
function,  Eq.  3.7,  into  the  Schrodinger  equation  gives 

H^ccsd  = Eccsd^ccsd-  (3-14) 

If  the  reference  energy  E0  is  subtracted  on  both  sides,  the  Schrodinger  equation  can 
alternatively  be  written  as 


Hn  ^ ccsd  — AEccsd'&ccsd,  (3.15) 

where  the  normal-ordered  Hamiltonian,  Eq.  3.6,  is  used.  A Eccsd  — ECcsd  — Eq 
is  the  correlation  energy  in  the  CCSD  model,  which  can  be  obtained  by  projecting 
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from  the  left  by  the  reference  function  <f>0 


AEccSD  = ($o|#Ar|'l> CCSd)- 


(3.16) 


The  intermediate  normalization  (Toi'b ccsd)  = 1 is  chosen.  Expanding  the  coupled- 
cluster  wave-function  leads  to 


A Eccsd  — ($o\Hn(1  + Ti  + T2  + -T^ITq). 


(3-17) 


Notice  that  the  exponential  series  is  self  terminating  because  of  the  two-particle 
character  of  the  Hamiltonian.  Recalling  that  the  vacuum  expectation  value  of  a 
normal-product  operator  is  zero,  the  energy  is  given  by 


^Eccsd  — ($o|#ivTi|<ho)  + ($o|-^vT2|$o)  E 2 
= L\  + L2  + L$. 


(3.18) 


{aUtfj}\$0)t%.  (3-19) 


As  an  example  the  term  L2  is  analyzed.  If  the  definitions  of  the  normal-ordered 
Hamiltonian,  Eq.  3.6,  and  the  T2  cluster  operator,  Eq.  3.9,  are  inserted,  L2  is  given 
by 

L2  = J2frsWs)  + ^^2{rs\\ut){rhhu} 

ijab  L rs  rstu 

The  contribution  of  the  one-electron  part  of  H ^ vanishes,  since  it  is  not  possible  to 
contract  all  the  operators  in  this  term  without  using  an  internal  contraction  in  one 
of  the  normal-ordered  products.  That  leaves 

L2  = ^^^(rs||ut)($o|{rtsttu}{ati6t;}|$0)t“6. 


ijab  rstu 


(3.20) 
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There  are  four  ways  to  contract  all  the  operators  between  the  normal  products 
£2  = 


ijab  rstu 


+ {rtV,*f°u00}{atV&t00.f  *}  + {rtV,,f°u00}{at0V*6t°.n 

T7T  ^ ^ ^ ^ (j"S 1 1 Ut')  (^ta ^ub^si ^rj  ^ta^ub^sj^ri  ^tb^ua^si^rj 


ijab  rstu 


+ fitbfiuafisjfiri)  ■ 


(3.21) 


In  addition  to  the  bullets,  circles  are  used  to  indicate  a certain  contraction.  Adding 
the  four  terms,  leads  to 

L2  = jI>'II  abK-  (3-22) 

ijab 

Similar  manipulations  can  be  done  for  the  L\  and  L3  terms,  which  gives  for  the 
CCSD  energy 

A Eccsd  = £ fia n + i £<y  ||o6)««»  + l £fo||ai)f ft*.  (3.23) 

ia  ijab  ijab 

The  coupled-cluster  energy  depends  on  the  amplitudes.  Equations  for  the  am- 
plitudes are  obtained  by  projecting  Eq.  3.15  onto  all  excited  determinants,  in  par- 
ticular onto  all  singly  and  doubly  excited  determinants  in  the  CCSD  model.  The 
number  of  resulting  equations  is  exactly  the  same  as  the  number  of  amplitudes  to 
be  determined.  The  amplitude  equations  for  the  CCSD  model  are 

A Eccsot*  = mnN  (l  + 7\  + r2  + h'l  + TjT2  + |4>o) 

a Eccsotj  = ($%\HN  [\  + Tl+  T2  + l-T{  + TXT%  + ^T22 

+ ^ + l-TlT,  + ^ |$0).  (3.24) 

An  important  feature  of  the  amplitude  equations  is  that  the  energy-dependent  terms 
on  the  left-hand  side  of  Eq.  3.24  cancel  with  disconnected  contributions  on  the 
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A 


V 


I II  III  W 

Figure  3-1.  One-electron  part  of  the  normal-ordered  Hamiltonian 

right-hand  side  (linked-diagram  theorem  [92,  93]).  The  amplitude  equations  can  be 
rewritten  as 

0 = W\HN(l  +T1+T2  + ±T?  + T1T2  + ±T?)\$o)c 

0 = (*$\HN(l  + T1+T2  + ±T*  + T1T2  + ±T* 

+ + \t?t2  + ±7?)  |$o)c,  (3.25) 

where  c stands  for  connected  contributions  only. 

It  is  possible  to  derive  the  amplitude  equations  using  Wick’s  thereom,  as  demon- 
strated above  in  the  example  of  the  L2  term  in  the  energy  expression.  However,  Eq. 
3.21  shows  that  those  manipulations  can  become  tedious.  For  that  reason  the  dia- 
grammatic formalism  is  introduced  next. 

The  diagrammatic  derivation  of  the  coupled-cluster  equations  is  a pictorial  way 
of  applying  Wick’s  theorem.  The  normal-ordered  Hamiltonian  and  the  cluster  op- 
erators in  diagrammatic  form  are  used  and  the  pieces  are  connected  appropriately, 
which  represents  a particular  class  of  contractions.  The  normal-ordered  Hamiltonian 
is  divided  into  a one-electron  and  a two-electron  part 

Hn  = Fn  + Wn-  (3.26) 

The  diagrammatic  representation  of  FN  is  given  in  Fig.  3-1,  and  its  algebraic 
interpretation  in  Table  3-1.  Down  going  lines  represent  holes,  up  going  lines  rep- 
resent particles.  Notice  that  in  the  Hartree-Fock  case  the  last  two  contributions  in 
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Table  3-1.  Algebraic  interpretation  of  FN 


diagram 

I II 

III  IV 

fabWb}  fij{Fj} 

faiWi}  fiai^a} 

A / 

\ 

A 

, + ' 

' + ' 

l 

A 7 

7\ 

\ 

> > 

> \ 

/ 

I II  III 


+ VV  - AA 

VIII  IX 


Figure  3-2.  Two-electron  part  of  the  normal-ordered  Hamiltonian 


Fig.  3-1  vanish  because  of  the  block-diagonal  nature  of  the  Fock  operator.  The 
diagrammatic  representation  of  W at  can  be  found  in  Fig.  3-2,  and  its  algebraic 
interpretation  in  Table  3-2. 

The  last  fragments,  needed  to  formulate  the  CCSD  equations  diagrammatically, 
are  the  amplitudes,  given  in  Fig.  3-3  and  in  Table  3-3.  In  order  to  obtain  the 
coupled-cluster  equations,  the  different  fragments  have  to  be  connected  in  all  possible 
ways  that  produce  a particular  excitation  level.  The  excitation  level  for  energy 
diagrams  is  zero,  for  7\  amplitudes  one,  and  for  T2  amplitudes  two. 

The  CCSD  energy  was  given  in  Eq.  3.18.  To  represent  the  first  term  in  Eq.  3.18, 
one  has  to  find  all  possible  ways  to  connect  a 7\  amplitude  with  the  Hamiltonian 
that  create  diagrams  with  zero  excitation  level  (i.e.  closed  diagrams,  where  all  lines 
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Table  3-2.  Algebraic  interpretation  of  W 


diagram 

I 

II 

III 

{ab\\cd){a)tfdc} 

(ij\\kl){i'jHk} 

(ai\\bj){aU^jb} 

diagram 

IV 

V 

VI 

VII 

(ab\\ci){a^bHc} 

(ia\\jk){i'tfkj} 

(ai\\bc)  {cdtf  cb} 

{ij\\ka){^fak} 

diagram 

VIII 

IX 

(ab\\ij){a)W  ji) 

(ij\\ab){i^fba} 

Tj  T2 

Figure  3-3.  Amplitudes 


lead  into  a vertex).  A vertex  is  the  Hamiltonian  or  an  amplitude.  There  is  only 
one  possibility  to  achieve  that.  The  second  term  has  to  be  represented  by  closed 
diagrams,  for  which  the  Hamiltonian  is  connected  to  a T2  amplitude.  Again  there 
is  only  one  way  to  do  that.  The  last  term  involves  closed  diagrams,  in  which  two 
Ti  amplitudes  are  connected  to  the  Hamiltonian.  Fig.  3-4  shows  the  CCSD  energy 
(algebraically  given  Eq.  3.23)  represented  in  diagrammatic  form. 

There  are  few  rules  to  follow,  when  evaluating  the  coupled-cluster  (antisym- 
metric) diagrams.  As  can  be  seen  in  Fig.  3-1,  Fig.  3-2,  and  Table  3-1,  Table  3-2, 
each  one-particle  Hamiltonian  vertex  is  associated  with  /out, in  > and  each  two-particle 
Hamiltonian  vertex  is  associated  with  (left  — out  right  — out || left  — in  right  — in). 
The  creation  and  annihilation  operators  are  implicitly  considered  by  connecting 
the  different  diagram  fragments  according  to  the  net  excitation  level  of  the  diagram. 
Internal  lines  are  summed  over,  those  are  lines  that  end  on  both  sides  in  a vertex. 
With  every  pair  of  equivalent  internal  lines  there  is  a factor  of  | associated.  A pair 
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Table  3-3.  Algebraic  interpretation  of  the  amplitudes 


diagram 

T,  T2 

of  internal  lines  is  equivalent,  if  they  connect  the  same  vertices  and  go  in  the  same 
direction.  There  is  also  a factor  of  | associated  with  every  pair  of  equivalent  T 
vertices.  T vertices  are  considered  equivalent  if  they  are  connected  with  the  same 
number  of  lines  equivalently  connected  to  the  Hamiltonian.  All  distinct  permuta- 
tions P of  inequivalent  external  lines  have  to  be  summed  over,  with  the  exception  of 
lines  that  would  become  equivalent  if  they  were  connected  by  a vertex.  External  lines 
are  lines  that  are  not  connected  on  both  sides  to  a vertex.  With  each  permutation 
a factor  of  (— where  cr(P)  is  the  number  of  necessary  exchanges  of  creation 
and  annihilation  operators  to  obtain  the  particular  permutation,  is  associated.  The 
sign  of  a diagram  is  determined  by  — where  h is  the  number  of  internal  hole 

lines  and  l is  the  number  of  loops. 

To  better  follow  which  diagram  belongs  to  which  term  in  the  amplitude  equa- 
tions, Eq.  3.25  is  divided  in  different  groups 

0 - (^l^^l  + Ti+Ta  + ^ + TxTa  + l^  |$0>c 

= S\  + S2  T <S3  + $4  + S5  + Sq 

o = ( i + r,  + t2  + iif  + TiT2  + ir| 

+ + ±T?T2  + ±T?)  |4„)e 

= Di  T D2  + D$  + D4  + + Dq  + Dy  + + Dg.  (3.27) 

Notice  that  the  connected  amplitude  equations  are  used,  where  disconnected  di- 
agrams were  already  canceled  by  energy-dependent  terms.  To  find  the  diagrams 
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Figure  3-4.  CCSD  energy 


contributing  to  the  equation  determining  the  7\  amplitudes  (all  S terms),  the  am- 
plitude vertices  have  to  be  connected  in  all  possible  ways  to  the  Hamiltonian  to  give 
a net  single  excitation,  where  a hole  line  and  a particle  line  is  not  connected  to  a 
vertex  from  the  top.  The  T2  amplitude  equations  can  be  obtained  by  finding  all 
possible  ways  to  connect  the  amplitudes  with  the  Hamiltonian  to  give  a net  double 
excitation,  where  two  hole  lines  and  two  particle  lines  are  not  connected  to  a vertex 
from  the  top.  In  Fig.  3-5  all  the  single  excitation  ( S ) terms  can  be  found.  The 
diagrams  for  the  T2  amplitudes  are  split,  in  Fig.  3-6  diagrams  Tfi  to  Z)4  are  given, 
in  Fig.  3-7  diagrams  D5,  and  in  Fig.  3-8  diagrams  D6  to  Dg.  The  diagrams  can 
be  interpreted  with  the  rules  given  above.  Notice  that  if  Hartree-Fock  orbitals  are 
used  diagrams  Si,  the  first  S3,  the  last  S4,  and  the  first  two  D5  are  zero. 

3.1.2  The  EOM-CC  Method 

Quantum  chemical  treatment  of  electron  correlation  within  the  coupled-cluster 
approximation  [90,  91]  traditionally  has  been  based  on  a single  Slater  determinant. 
That  usually  limits  its  applicability  to  the  lowest  electronic  states  of  a particu- 
lar symmetry  and  multiplicity.  An  efficient  tool  to  access  electronically  excited 
states  is  the  equation-of-motion  coupled-cluster  (EOM-CC)  method  [94-96].  The 
excited-state  wave-function  is  generated  from  the  coupled-cluster  wave-function  by 
the  action  of  a wave  operator.  The  wave  operator  is  not  restricted  to  conserve  the 
number  of  particles.  If  the  excited  state  is  chosen  to  be  a state,  where  an  electron 
has  been  removed,  the  IP-EOM-CC  method  [97]  emerges.  If  an  electron  is  added, 
while  forming  the  excited  state,  the  approach  is  referred  to  as  EA-EOM-CC  method 


53 


Figure  3-5.  Single-excitation  CCSD  equations 


[98].  Here,  the  focus  is  on  the  IP-EOM-CC  and  the  EA-EOM-CC  methods,  since 
those  are  reformulated  to  obtain  a correlation  potential  for  the  CIP  model. 

The  equation  of  motion  is  given  by 

[Hn,  fJfc]  = ^k^k^o-  (3.28) 

is  the  reference  wave- function  and  HN  is  the  normal-ordered  operator,  introduced 
before  in  Eq.  3.6.  The  operator  ttk  consist  of  a constant  and  a string  of  creation  and 
annihilation  operators.  If  it  preserves  the  number  of  particles,  a >k  is  an  excitation 
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D4  Dj  D4 

Figure  3-6.  Double-excitation  CCSD  equations,  terms  D\  through  D,\ 

energy  Ek  — E0  = A Ek  — AE0  = ujk.  However,  flk  may  as  well  change  the  number 
of  particles,  making  different  sectors  of  the  Fock  space  accessible. 

The  ( N — l)-particle  eigenstates  may  be  obtained  by  using 

ttk  = YlCi  W + cij  {^5*}  + • • • • (3-29) 

i 

In  this  case  tok  is  an  ionization  energy.  The  (N  + l)-particle  eigenstates  can  be 
accessed  by  using 

nk  = J2 °a  {at}  + E °ib  ft™'}  + ■■■■  (3.30) 

a i,a<b 

Choosing  the  wave  operator  as  in  Eq.  3.30,  causes  uk  to  be  the  negative  of  an 
electron  affinity.  The  braces  refer  to  normal  ordering.  The  parameters  c\  are  the 
amplitudes  to  be  determined  for  the  linear  expansion  of  the  eigenstates. 

In  the  EOM-CC  approach  the  reference  function  is  chosen  to  be  the  coupled- 
cluster  wave- function,  Eq.  3.7,  specifically,  in  the  EOM-CCSD  model,  the  CCSD 
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D5  d5 

Figure  3-7.  Double-excitation  CCSD  equations,  terms  D5 

wave- function  3.13  is  used.  Usually  a balanced  description  for  the  reference  function 
and  the  linear  expansion  of  eigenstates  is  employed,  which  means  that  in  the  IP- 
EOM-CCSD  and  the  EA-EOM-CCSD  model  at  most  a single  excitation  in  addition 
to  an  ionization  or  attachment  process  are  considered  (i.e.,  the  dots  in  Eq.  3.29  and 
Eq.  3.30  disappear). 

Inserting  ansatz  3.7  in  Eq.  3.28,  the  equation-of- motion  becomes 

[Hn,  fife]  eT<ho  = <^fcDfceT<f>0.  (3.31) 

Multiplying  on  the  left  with  e-T,  and  and  using  the  fact  that  [f2,  T]  — 0,  leads  to 

[Hn,  fyfc]  $o  = (3.32) 

where  HN  is  the  transformed  Hamiltonian,  consisting  only  of  the  connected  part  of 
the  operator  product 


Hn  = e~THNeT  = ( HNeT)c . 


(3.33) 
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Hn  contains  all  the  ground-state  CC  information. 

Eq.  3.32  can  be  rewritten  into  an  eigenvalue  problem  by  considering  the  fol- 
lowing 

[Hni  fifc]  $o  = ET/v^fcd>0  — 

= — f^fcAE'c'c^o 

= (AEfc  — AEccO^fc^o 

= (Hn  — A£'c'c-)f2fc<f>o-  (3.34) 

Since  A Ecc  represents  a closed  diagram,  Eq.  3.34  can  also  be  written  as 


[Hn,  fife]  d*0  = (HN^k)c  $0- 


(3.35) 
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Therefore,  Eq.  3.32  can  be  formulated  as  a nonhermitian  Cl-like  eigenvalue  problem 
for  each  eigenvalue  uk. 

Because  of  the  nonhermiticity  of  Hn,  there  is  a right  and  a left  eigenvalue 
problem,  where  the  left  and  right  eigenvectors  correspond  to  the  same  eigenvalues. 
In  matrix  form  Eq.  3.32  is  given  by 

(HjvCfc)c  = Ckcuk,  C[,Hjv  = C[ujk.  (3.36) 


The  eigenvectors  Ck  contain  the  amplitudes  C\  from  Eq.  3.29  or  Eq.  3.30,  depending 
on  whether  Hat  is  expressed  in  the  basis  of  hole  determinants,  corresponding  to  the 
(N  — l)-particle  states,  or  in  the  basis  of  particle  determinants,  corresponding  to 
the  (./V  + l)-particle  states.  The  eigenvalues  uk  are  the  exact  ionization  potentials  or 
the  negative  of  the  exact  electron  affinities  respectively,  if  no  truncation  is  applied. 

If  the  excitation  space  is  chosen  according  to  the  IP-EOM-CCSD  model  or  the 
EA-EOM-CCSD  model,  the  ionization  potentials  or  electron  affinities  are  in  general 
well  approximated.  Hat  in  Eq.  3.36  then  takes  the  form 


H N = 


( (S\Hn\S)  (S\Hn\D)  ^ 

K ( d\hn\s ) ( d\hn\d ) J 


(3.37) 


where  S stands  for  determinants  with  a single  ionization  or  attachment  process,  and 
D for  determinants  with  an  additional  excitation. 

The  left  and  right  eigenvectors  of  a nonhermitian  matrix  form  a biorthogonal 
set,  which  can  be  chosen  to  be  biorthonormal,  so  that 


C(A)  • C(m)  = clX)  • = V (3.38) 

V 

For  simplicity  the  subscript  k,  discriminating  between  the  different  choices  of  ff  is 


omitted. 


58 


Figure  3-9.  H/p  for  ionization  potentials 


Since  the  wave-function  is  explicitly  known  in  the  EOM-CC  method 

<*aI  = £4A)<4>,|  l*„>  = (3.39) 

V K 

properties  can  be  calculated  as  a generalized  expectation  value 

6 = Tr  (DNQ),  (3.40) 

where  © is  an  arbitrary  linear  operator.  The  N-particle  reduced-density  matrix  is 
defined  as 

Dpq... rs  = (®a|pV  • - • sr\yx).  (3.41) 

However,  the  N-particle  reduced-density  matrix  is  not  size  consistent  [99],  and  as  a 
result,  properties  calculated  as  expectation  values  depend  on  the  physical  extent  of 
the  system.  Properties  should  preferentially  be  defined  as  energy  derivatives. 
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p Q 


Figure  3-10.  \HEA  for  electron  attachment 


The  particular  realization  of  H for  the  IP-EOM-CCSD  model  is  diagrammat- 
ically  shown  in  Fig.  3-9,  and  for  the  EA-EOM-CCSD  model  in  Fig.  3-10.  The 
space  is  conveniently  partitioned  into  the  P space,  which  includes  in  the  case  of 
IP-EOM-CCSD  all  determinants  resulting  from  a single  excitation  of  an  electron 
into  the  continuum,  also  called  l-hole  space,  and  in  the  case  EA-EOM-CCSD  all 
determinants,  resulting  from  the  attachment  of  one  electron,  also  called  1-particle 
space.  In  the  IP-EOM-CCSD  model  the  Q space  consists  of  all  determinants,  result- 
ing from  an  ionization  in  addition  to  a single  excitation,  the  2-hole  1-particle  space. 
In  the  case  of  EA-EOM-CCSD  the  Q space  is  composed  out  of  all  determinants, 
obtained  by  an  attachment  in  addition  to  a single  excitation,  the  2-particle  1-hole 
space.  Wiggling  lines  in  Fig.  3-9  and  Fig.  3-10  refer  to  H interactions.  Their  in- 
terpretation can  be  found  in  Fig.  3-11.  Diagonalization  of  H/F  delivers  the  normal 
IP’s  in  the  PP  block,  and  the  shake-up  IP’s  in  the  QQ  block.  Diagonalization  of 
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F igure  3-11.  Interpretation  of  H interactions 


YLEA  yields  the  negative  of  the  normal  EA’s  in  the  PP  block,  and  the  negative  of 
the  shake-up  EA’s  in  the  QQ  block. 

3.1.3  Partitioned  EOM-CC  Approach 

In  this  section  |a  correlation  potential  that  yields  the  exact  IP’s  and  the  exact 
EA’s  as  the  negative  of  the  orbital  energies  is  discussed.  The  EOM-CC  method  is 
exploited,  since  it  gives  the  exact  IP’s  and  EA’s  if  not  truncated. 
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Figure  3-11.  Continued 


H/p  and  HEA  as  given  in  Fig.  3-9  and  Fig.  3.10  cannot  directly  be  used  to 
extract  a correlation  potential,  since  the  operators  act  in  their  P+Q  spaces.  To 
obtain  a correlation  potential  for  an  independent  particle  model  according  to  Eq. 
2.96,  the  operators  have  to  act  only  in  their  P spaces.  Matrix  partitioning  [100] 
can  be  applied  to  formulate  an  effective  Hamiltonian  with  reduced  dimension.  This 


Figure  3-12.  Illustration  of  Vc 


well-known  technique  has  been  used  to  simplify  excited  state  calculations  [101,  102], 
and  calculations  of  NMR  coupling  constants  [103]. 

The  partitioning  of  H/p  and  YlEA  in  P and  Q space  was  introduced  in  Fig.  3-9 
and  Fig.  3-10.  The  eigenvectors  in  Eq.  3.36  can  be  partitioned  accordingly,  leading 
to  the  partitioned  equation  of  motion.  However,  since  the  IP-EOM-CC  method 
yields  IP’s,  but  the  orbital  energies  of  the  CIP  model  are  sought  to  be  the  negative 
of  the  IP’s,  the  partitioned  equation  is  multiplied  by  minus  one.  The  negative  of  the 
IP  is  denoted  by  e.  The  partitioned  IP-EOM-CC  equation  is  then  given  by 


(3.42) 


Eq.  3.42  can  be  divided  into  two  equations 


(3.43) 


(3.44) 


Solving  Eq.  3.44  for  C q leads  to 


(el  + H"^-1  H'TpCV  = Ce. 


(3.45) 


Substituting  Eq.  3.45  into  Eq.  3.43  gives 


-H£  + H^(el  + H^)_1H'£  CP  = £Cp. 


(3.46) 
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An  effective  Hamiltonian  H can  now  be  defined 

H = — Hpp  + Hp^  (el  + H&)  H‘qpp,  (3.47) 

which  yields  the  negative  of  the  principal  IP’s  upon  diagonalization.  A new  eigen- 
value equation  is  obtained,  which  is  solved  only  in  the  P space,  where  the  effects  of 
the  Q space  are  included  through  the  effective  Hamiltonian 

&eff(£)Cp  = sCp.  (3.48) 

Partitioning  YiEA,  Fig.  3-10,  in  the  basis  of  net  particle  determinants  leads 
to  the  effective  Hamiltonian  that  yields  the  negative  of  the  principal  EA’s  upon 
diagonalization 

Hf/,(e)  = Hpp  + HpQ  (el  - fl§3) H“.  (3.49) 

There  is  no  need  to  multiply  by  minus,  like  for  the  IP-EOM  model.  The  effective 
eigenvalue  problem  solved  in  the  1-particle  space,  is  given  by 

Hf//(e)CP  = eCP.  (3.50) 

Notice  that  the  Q spaces  are  not  restricted  to  the  2- hole  1-particle  space  or 
the  2-particle  1-hole  space  respectively,  but  can  consist  of  any  number  of  particle 
and  hole  combinations  that  amounts  to  a net  hole  determinant  or  a net  particle 
determinant.  If  no  restrictions  to  the  Q space  are  applied,  the  eigenvalues  of  the 
effective  Hamiltonians  are  the  negative  of  the  formally  exact  principal  IP’s  or  the 
negative  of  the  formally  exact  principal  EA’s. 

If  the  effective  Hamiltonians  are  expressed  in  canonical  Hartree-Fock  orbitals, 
then  the  diagonal  elements  contain  the  Hartree-Fock  orbital  energies.  This  can 
be  seen  by  expanding  the  diagrams  in  the  PP  blocks  according  to  Fig.  3-11.  The 
leading  contribution  to  [H iff]ij  is  /tJ  = in  canonical  Hartree-Fock  orbitals.  The 
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leading  contribution  to  [H ffyab  is  fab  = ea(Ab  in  canonical  Hartree-Fock  orbitals. 
The  diagonal  elements  of  the  effective  Hamiltonians  contain  a correlation  correction 
as  well,  the  off-diagonal  elements  contribute  only  to  correlation.  Eq.  3.48  and  Eq. 
3.50  can  be  written  as 

H£(e)C">  = (eIP  + VIcp(e))C'Pp  = eCIPp  (3.51) 

HffAf(e)CpA  = (eEA  + VEA{e))CfA  = £Cf,A,  (3.52) 

with 

Vf(e)  = — (HppJcOTr  + Hpp  (el  + Hqq)~1  (3.53) 

Vf*(e)  = + (3.54) 

Eq.  3.53  and  Eq.  3.54  define  a correlation  potential  in  the  1-hole  and  the  1-particle 
space.  To  obtain  a correlation  potential  in  the  unified  P spaces,  'VIcp(e)  and  YEA(e) 
are  added  as  shown  in  Fig.  3-12.  The  eigenvalue  equation  in  the  CIP  model,  derived 
from  the  partitioned  EOM-CC  method,  is  given  by 

He//(e)C(e)  = (e  + Vc(e))C(e)  = eC(e).  (3.55) 

Notice,  even  though  Eq.  3.55  yields  the  negative  of  the  exact  IP’s  and  EA’s  if  no 
truncation  is  introduced,  there  is  a certain  arbitrariness  if  one  is  interested  in  the 
eigenvectors.  Since  the  IP-EOM-CC  and  the  EA-EOM-CC  methods  are  completely 
independent,  the  off-diagonal  blocks  of  the  correlation  potential  are  not  defined  and 
set  to  zero.  However,  if  one  of  the  off-diagonal  blocks  were  not  equal  to  zero,  the 
same  IP’s  and  EA’s  would  be  obtained,  but  the  eigenvectors  would  be  different. 
That  point  is  discussed  further  in  Section  3.3. 

Notice  also,  in  Eq.  3.55  e appears  not  only  as  the  eigenvalue,  but  Vc(e)  also 
depends  on  e,  which  in  turn  causes  the  eigenvectors  to  depend  on  e.  That  means 
that  for  each  IP  or  EA  a new  eigenvalue  problem  has  to  be  solved. 
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The  effective  Hamiltonians  in  Eq.  3.51  and  Eq.  3.52  act  in  the  space  of  N- 
electron  Hartree-Fock  determinants,  from  which  an  electron  has  been  removed,  if 
the  IP  sector  is  considered,  or  to  which  an  electron  is  added,  if  the  EA  sector  is 
considered.  However,  each  element  of  H ^ is  defined  by 

= (<h\R‘4M>s),  (3.56) 

where  the  \(j>i)  are  one-electron  functions,  represented  by  C pp,  solved  for  the  IP 
sector.  Equivalently,  each  element  of  H fA  is  defined  by 

[Hf^U  = ($„|a5“6'|4.„>  = (3.57) 

Again,  the  \<f>a)  are  one-electron  functions,  represented  by  C pp,  solved  for  the  EA 
sector. 

He//  = e + Vc  can  equivalently  be  regarded  as  an  effective  one-particle  operator 
in  an  orbital  space,  where  the  orbitals  are  solutions  to  the  eigenvalue  equation 

= (/  T Uc(£m)) |0m)  = £m|0m) • (3.58) 

/ is  the  Fock  operator,  and  vc  is  the  correlation  potential  given  in  matrix  form  in 

Fig.  3-12. 

The  effective  Hamiltonian  for  the  EOM-CCSD  model  was  implemented.  Eq. 
3.55  can  be  solved  iteratively,  starting  with  Koopmans’  IP’s  and  EA’s  as  a guess. 
The  procedure  was  tested  on  various  small  molecules,  and  convergence  to  the  IP- 
EOM-CCSD  and  the  EA-EOM-CCSD  values  was  easily  achieved  in  most  cases. 

The  possibility  of  approximating  the  effective  Hamiltonians  is  also  considered. 
The  simplest  approximation  is  He//  = Hpp  and  Hfyj  = Hfp  respectively.  Test 
calculations  on  simple  systems  showed,  that  in  most  cases  this  approximation  is  not 
even  an  improvement  over  Hartree-Fock.  As  one  would  expect,  the  influence  of  the 
Q space  in  the  effective  Hamiltonian  is  of  great  importance  and  cannot  be  neglected. 
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Table  3-4.  Comparison  of  different  approximation  schemes  for  IP’s  and  EA’s 
for  H2,  He2  and  HF,  IP’s  and  EA’s  are  given  in  Hartree 


HF 

I 

II 

III 

IV 

H2  IP  (DZP) 

0.59299 

0.62792 

0.60167 

0.60157 

0.60137 

H2  EA  (DZP) 

-0.26004 

-0.83599 

-1.17694 

-0.26912 

-0.84292 

-1.18550 

-0.24775 

-0.81684 

-1.10297 

-0.24872 

-0.81452 

-1.02271 

-0.24885 

-0.81595 

-1.09046 

He2  IP  (DZ) 

0.91595 

0.91231 

0.93153 

0.92805 

0.88355 

0.87986 

0.87738 

0.87354 

0.87845 

0.87462 

He2  EA  (DZ) 

-1.37551 

-1.42489 

-1.39015 

-1.43918 

-1.35772 

-1.40595 

-1.35060 

-1.39845 

-1.35118 

-1.39910 

HF  IP  (DZ) 

0.64285 

0.75240 

0.68434 

0.80442 

0.53925 

0.68949 

0.55274 

0.70360 

0.55843 

0.70596 

HF  EA  (DZ) 

-0.21248 

-1.05750 

-1.12680 

-0.22614 

-1.07967 

-1.14724 

-0.20102 

-1.01360 

-1.03460 

-0.20265 

-1.02724 

-1.03070 

-0.20278 

-1.02838 

-1.04402 

Some  illustrative  results  are  summarized  in  Table  3-4.  The  column  denoted  by  I 
contains  results  of  the  approximate  scheme  just  mentioned,  the  Hartree-Fock  results 
are  included  for  comparison  in  the  column  denoted  by  HF,  and  the  last  column, 
denoted  by  IV,  shows  results  of  the  fully  iterated  Eq.  3.55,  which  coincide  with  the 
EOM-CCSD  results  for  IP’s  and  EA’s. 

Since  the  Q space  is  of  much  higher  dimension  than  the  P space,  it  would  be 
useful  to  be  able  to  simplify  Hqq  an<^  Hqq-  Since  the  diagonal  elements  of  H(/g 
and  Hqq  are  larger  than  the  off-diagonal  elements,  the  approximation  of  H and 
Hqq  by  only  its  diagonal  elements  was  tested.  Results  for  the  test  examples  are 
included  in  Table  3-4  in  the  column  denoted  by  II.  The  results  are  good. 

It  was  also  tested  if  the  iteration  of  Eq.  3.55  is  necessary,  or  if  good  results 
can  be  achieved  by  using  Koopman’s  IP’s  and  EA’s  in  Eq.  3.53  and  Eq.  3.54  for  e 
without  iterating.  Results  of  this  approximation  are  included  in  Table  3 4 in  column 
III.  The  test  calculations  showed  that  the  approximations  II  and  III  for  Hqq  and 
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H qq  are  quite  reasonable  and  can  be  considered  as  possible  candidates  to  extract 
correlation  potentials  from  them. 

With  the  exception  of  the  simplest  approximation,  which  is  not  recommended, 
the  above  schemes  all  have  in  common  that  the  correlation  potential,  resulting  from 
them,  depends  on  e.  That  means  that  for  each  orbital  a different  correlation  poten- 
tial would  have  to  be  used,  and  that  each  IP  and  EA  has  to  be  calculated  separately. 
It  also  implies  that  the  orbitals  are  not  orthogonal.  For  those  reasons  the  energy 
dependence  of  the  correlation  potential  is  viewed  as  a disadvantage,  and  in  a later 
section  an  energy-independent  correlation  potential,  which  is  universal  for  all  or- 
bitals, is  discussed. 

3.2  Electron-Propagator  Method 

An  alternative,  formally  exact  one-electron  theory,  can  be  derived  from  the 
electron-propagator  [18].  In  contrast  to  the  EOM-CC  methods,  the  1 hole  and  the 
1 particle  sectors  are  not  decoupled,  but  simultaneously  described.  The  self-energy 
plays  the  role  of  the  correlation  potential,  which  incorporates  all  the  many-body 
effects.  The  self-energy  is  generally  energy  dependent,  which  causes  the  orbitals  to 
be  energy-dependent  as  well.  They  are  therefore  linearly  dependent. 

In  the  treatment  of  electron  propagators,  it  is  common  to  use  the  so-called 
Heisenberg  picture  [89].  To  distinguish  the  wave- function  and  the  operators  from 
the  usual  Schrodinger  picture,  the  superscript  S is  used  for  the  Schrodinger  represen- 
tation, and  H for  the  Heisenberg  representation.  The  time-dependent  Schrodinger 
equation  becomes 

ih^s{t))^Hs\<Hs{t)),  (3.59) 

where  Hs  is  assumed  to  have  no  explicit  time  dependence.  Since  Eq.  3.59  is  a first- 
order  differential  equation,  the  initial  state  t0  determines  the  subsequent  behavior. 
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A formal  solution  to  Eq.  3.59  is  given  by 

\yS(t))  = ez*^\*s(t0)). 


(3.60) 


Since  H is  hermitian,  the  exponential  represents  a unitary  operator 


Tr.  . -iHs(t-tn) 

U(t)  = e » 


(3.61) 


called  the  evolution  operator,  which  satisfies 

lft(t)U(t)  = U(t)U\t)  = 1.  (3.62) 

The  evolution  operator  carries  the  initial  state  |^,5(to))  into  the  final  state  |'I's(f)). 
If  \Vs(t0))  is  a time-independent  eigenfunction  of  Hs , the  exponential  in  Eq.  3.60 

—iEt 

becomes  the  factor  e~~ for  a stationary  state  of  energy  E. 

To  obtain  the  wave-function  and  operators  in  the  Heisenberg  picture,  a unitary 
transformation  is  applied  according  to 

|*H)  = U*(t)\Vs{t))  AH(t)  = U*(t)As{t)U(t).  (3.63) 


The  wave-function  in  the  Heisenberg  representation  is  time  independent,  since 


I*")  = u'(t)mt)\<s>s{t0))  = |'KS((„)). 


(3,64) 


Operators  in  general  acquire  a time  dependence.  Substituting  the  definition  of  the 
evolution  operator  into  Eq.  3.63,  a general  Heisenberg  operator  is  given  by 


AH(t)  = e"  » As{t)e 


iHS  (t—tft)  q - iHS(t-tn ) 

t A ^ ! 4-  \ ^ t 


(3.65) 


If  A5  is  time  independent,  the  time  derivative  of  Eq.  3.65  yields 


t ^ iH(t-tn)  q —iH(t—to)  jj 

* [As,H]e  * =[AH(t),H ], 


(3.66) 
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where  the  superscript  in  the  Hamiltonian  is  omitted,  since  H means  the  usual  time- 
independent  Hamiltonian.  Eq.  3.66  establishes  the  equation-of-motion  of  an  opera- 
tor in  the  Heisenberg  picture.  In  particular,  if  As  commutes  with  H , the  right-hand 
side  vanishes,  and  AH  is  a constant  of  the  motion. 

The  Green’s  function,  also  called  propagator,  is  an  expectation  value  of  field 
operators.  The  N-particle  many-body  Green’s  function  is  defined  in  the  Heisenberg 
representation  as 

Gpq...sr(tp,tq  • ■ -ts,tr)  = (— i)N \T\pH  {tp)r^H  (tr)qH (tq)s^H(ts) . . .] (3.67) 


where  the  operators  pH(tp ) and  A11  (tr)  are  the  annihilation  and  creation  operators 
in  the  Heisenberg  picture,  with  the  time  dependence  given  by 


vH(tP) 


iH(t-t  n)  -iH(t-tn) 

— e » p e a 


rW 


(tr)  — 


iH(t-tn)  . -iH(t-tn) 

e * r'e  h 


(3.68) 


\^q)  is  the  Heisenberg  ground  state  of  the  interacting  system,  which  is  assumed 
to  be  normalized  = 1 . The  T product  of  several  operators  orders  them 

from  right  to  left  in  ascending  time  order  and  adds  a factor  of  (— l)p,  where  P is 
the  number  of  interchanges  of  fermion  operators  from  the  original  given  order.  The 
N-particle  Green’s  function  contains  N pairs  of  creation  and  annihilation  operators, 
each  having  a unique  time  argument. 

The  single  particle  Green’s  function  contains  only  one  pair  of  creation  and 
annihilation  operators 


(3.69) 
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Taking  the  T ordering  explicitly  into  account,  Eq.  3.69  can  be  written  as 


W\pH(tPWH(tq)W)  tp  > tq 

{^oWHitq)pH{tp)\^o)  ^ > tp 


(3.70) 


Eq.  3.70  can  be  conveniently  expressed  by  using  the  Heaviside  step  function  9 


These  states  are  eigenstates  of  the  Hamiltonian  and  include  all  possible  numbers  of 
particles 


The  superscript  H is  omitted,  since  the  stationary  solutions  to  the  Schrodinger  equa- 
tion satisfy  the  time-independent  form  of  the  Schrodinger  equation.  The  integrals 
in  Eq.  3.73  can  only  be  nonzero  if  the  states  |^n)  contain  N ± 1 particles,  while  the 
state  I'ko)  contains  N particles. 

The  Green’s  functions  are  often  studied  in  terms  of  their  Fourier  transforms 


Gpq(tp,tq)  = -i9(tp-tq)(<b”\pH(tp)q]H(tq)\^) 

+ i9(tq  - tp){^\q^H (tq)pH (tp)\^H). 


(3.71) 


A complete  set  of  Heisenberg  states  is  inserted  between  the  field  operators. 


n 


+ i»((,-g«i<;»''(t,)i'p")('pnH|p"(giO].  (3.72) 


Substituting  Eq.  3.68  into  3.72  gives 


Gpq(tp,  tq ) — ^ ^ I"  i9(tp  tq) 


— i(En  — Eo)(tp  — tg) 

e * 


(^o|p|^)(^|gt|^o) 


n 


— i(EQ  — En)(tp  — tq) 

h 


('J'ol«,|4,„)(4,„|p|'I'o)  • (3.73) 


(3.74) 
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with  the  inverse  relation 


1 r°° 

G(t,t')  = — G{E)e-iE(t~tf)dE. 
2tt 


(3.75) 


Using  the  identities 


hm  -5-  [ 

v-*o  2n  J _ 


1 r°°  e-iE(t-t') 


-dE  = 


t > t' 


E-(E0-  En ) - IT]  j ie-i(E0-En)(t-t')  t <ti 


(3.76) 


and 


lim—  f 
v-o  2n  J_c 


-iE(t-f)  I -ie-i(En-Eo)(t-t>)  t > t, 

-dE  = < (3.77) 


E - ( En  - E0)  + irj 


0 


t < t' 


The  Fourier  transform  can  be  written  as 


Gpq(E)  = Ihri  ^ 


(^0|plT»)(Tn|gt|^0)  ^ (^0|gt|^w)(^n|p|^0)- 
E — En  + Eq  + if]  E — Eq  + En  — irj 


(3.78) 


This  is  the  so  called  spectral  or  Lehmann  representation  of  the  propagator.  A com- 
mon notation  for  Gpq(E)  is  also  ((p;  q^))E-  The  electron-propagator  in  its  spectral 
representation  is  energy  dependent,  and  poles  occur  when  the  energy  E is  equal  to 
the  negative  of  an  electron  affinity  -( E0(N ) - En(N  + 1))  or  the  negative  of  an 
ionization  potential  —(En(N—  1)  — E0(N)). 

The  starting  point  for  most  approximation  schemes  is  the  equation-of-motion, 
which  can  be  solved  directly,  without  the  need  to  seek  a variational  approximation 
to  the  many-electron  wave-function.  Using  the  identity 


E(E  — a)  1 = \ + a(E  — a) 


(3.79) 


and  the  relations 


('I'o|p|'I'n)(£n  - £o) 

- Bo) 


(3.80) 
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two  equivalent  forms  of  the  equation-of- motion  can  be  derived  from  Eq.  3.78 

E({p\q]))E  = ('J'oltzWl+l^o)  + ((\p,H]-,q*))E 

= (^olb.^l+l^o)  + {(p;[H,qi]))E.  (3.81) 

The  so  called  moment  expansion  can  be  obtained  by  iterating  the  equation-of- 
motion.  Iterating  the  second  of  Eq.  3.81  gives 

((?;?*»£  = B'1('I'ol[p.9,]+|4'o>  + £;‘2('I'ol[p.[ff,9,]]+l'I'o) 

+ £r»(*o||p,  \H,  [H,  <jt]]]+|9'„)  + . . . . (3.82) 

To  represent  the  propagator  compactly,  and  to  use  standard  matrix  and  vector 
space  techniques,  superoperators  can  be  employed  [104,  105].  A linear  space  is 
introduced,  whose  elements  are  linear  combinations  of  field  operator  products 

{z,  a}  U {aUj,  i^ab}  U (3.83) 

Superoperators  act  on  this  space,  and  the  superoperator  Hamiltonian  and  the  su- 
peroperator identity  are  defined  by 

HX  = [X,  H] 

IX  = X,  (3.84) 

where  A is  a general  element  of  the  linear  space.  The  scalar  product  (X|y),  where 
Y is  another  element  of  the  linear  space,  is  defined  as 

(X|y)  = (^o|[y,A-t]+|$0).  (3.85) 

The  moment  expansion,  Eq.  3.82,  can  then  be  expressed  as 


(( V\q]))E  = E \q\p)  + E 2(q\Hp)  + E 3{q\H2p)  + . . . . 


(3.86) 
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A formal  summation  is  performed  to  obtain 

((P'<Q*))e  = (3.87) 

where  ( EI—H)~l  is  the  superoperator  resolvent.  Introducing  vector  arrays  of  simple 
field  operators  a,  the  electron-propagator,  Eq.  3.87,  can  be  written  in  matrix  form 

G(E)  = (a\(Ei  - H)-1*).  (3.88) 

The  superoperator  resolvent  (El  - H)~l  can  be  reformulated  as  a matrix  in- 
verse by  inner  projection  [106].  With  an  orthonormal  truncated  basis  e,  there  is  a 
projection  operator  p associated 

P = le*Xeil  = ee+-  (3.89) 

i 

If  the  basis  is  not  orthonormal,  the  metric  matrix  S = e^e  must  be  included.  Using 
instead  of  e the  Lowdin  basis  e = eS  2 , it  follows 

P = le»)(‘5'-1)ii(ejl  = eS_1ef.  (3.90) 

ij 

For  a complete  space,  p becomes  the  unit  operator.  When  a truncated  basis  is  used 
in  operator  space,  there  are  two  kinds  of  useful  projections,  the  outer  and  the  inner 
projection 

A!  = pAp  A"  = A^pA%.  (3.91) 

The  inner  projection  is  particularly  useful  in  propagator  theory.  Neither  projec- 
tion preserves  the  relationship  AB  = C.  Nevertheless,  the  projected  forms  may 
be  regarded  as  approximations,  which  become  increasingly  precise,  as  the  basis  is 
extended. 

Inserting  Eq.  3.90  in  the  inner  projection  form  of  Eq.  3.91,  leads  to 


A"  = AseS'VAU 


(3.92) 
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A is  assumed  to  be  Hermitian  and  positive  definite.  Introducing  the  transformed 
basis  u 

A^e  = (A^e\  A^e2  . . .)  = u,  (3.93) 

the  metric  becomes 

S = e*e  = = uM_1u.  (3.94) 


Using  Eq.  3.93  and  Eq.  3.94,  the  inner  projection  Eq.  3.92  can  be  written  as 

A"  = u(utA-1u)-1ut.  (3.95) 

Since  A is  an  arbitrary  operator,  it  may  be  replaced  by  A-1.  Defining  ifiAu  = A, 
with  Aij  = (vn\A\uj),  the  inverse  of  the  operator  A can  be  approximated  in  terms 
of  a matrix  inverse 

(A"1)"  = u(utAu)-1ut  = uA-V.  (3.96) 

If  A is  taken  to  be  the  superoperator  resolvent  (El  — H)~l,  and  u is  a vector 
of  all  electron  field  operators,  Eq.  3.88  can  be  rewritten  as 

G(E)  = (a|u)(u|(£7  - jy)u)_1(u|a),  (3.97) 


where  the  inverse  of  a matrix  instead  of  the  resolvent  operator  is  employed.  If  the 
manifold  u spans  the  entire  space,  no  approximation  is  made. 

u can  be  partitioned  into  the  primary  space  of  simple  field  operators  a and  the 
orthogonal  complement  f.  Eq.  3.97  then  becomes 


G (E)  = 


1 0 


/ „ \ 

El  - (a|77a)  -(a|77f) 

^ -(f|^a)  El  - (f|£f)  } 


< 1 ' 


\ 


(3.98) 


where,  because  of  the  orthogonality  condition  (a|f)  = 0,  only  the  upper  left  corner 
of  the  inverse  matrix  contributes. 
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The  manipulations  to  obtain  the  inverse  of  Eq.  3.98  are  done  symbolically  for 
simplicity.  Upper  case  letters  are  used  for  the  elements  of  the  matrix  in  Eq.  3.98, 
and  lower  case  letters  are  used  for  the  inverse.  The  condition  a matrix  and  its  inverse 
have  to  fulfill  is 


/ 


A B 


\ 


\ 


'l  0^ 


(3.99) 


V 


C D 


\ 


c d 


\ 


0 1 


From.  Eq.  3.99  the  following  equations  can  be  obtained 


Aa  + Bc  = 1 (3.100) 

Ca  + Dc  = 0.  (3.101) 

Solving  Eq.  3.101  for  c leads  to 

-D-'Ca  = c,  (3.102) 

which  is  substituted  in  Eq.  3.100  and  solved  for  a 

a=  (A-  BD~lC)~l . (3.103) 

Substituting  the  appropriate  terms  from  Eq.  3.98  into  Eq.  3.103,  yields  the  parti- 
tioned form  of  the  inverse  propagator  matrix 

G-1(jE)  = (a|(£7  - H) a)  - (a|tff)(f|(£J  - ^fj-^f^a).  (3.104) 


Various  approximate  schemes  can  be  introduced  by  truncating  the  manifold  f 
and  by  different  choices  of  the  approximate  ground  state  |T0),  which  defines  the 
scalar  product,  Eq.  3.85.  H is  commonly  chosen  to  be  separated  in  the  unperturbed 
superoperator  H0  and  the  perturbation  superoperator  V 

ho  = Y1  ep(PjP}- 

v 


H = H0  + V, 


(3.105) 
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With  the  definition 

Go1(£)  = (a|(E/-ff0)a)>  (3.106) 

the  Dyson  equation  can  be  written  as 

G~1(E)  = Go1{E)-'E(E).  (3.107) 

The  self-energy  'E(E)  is  that  part  of  G-1(.E)  that  contains  dynamical  relaxation 
and  correlation  effects 

£(£)  = (a|Va)  + (a\Hf)[El  - (f | ^f)]-1  (f |^a).  (3.108) 

As  already  mentioned,  IP’s  and  EA’s  are  given  by  the  negative  of  the  poles  of 
G(E).  When  detG(E')  — ► oo,  or  equivalently,  when  detG_1(£l)  — * 0,  E is  a pole.  It 
is  equivalent  to  requiring  G~1(E)  to  have  a vanishing  eigenvalue  at  the  pole  energy 

G~l(E)C(E)  = 0 C(E).  (3.109) 

If  |tf0)  is  defined  to  be  the  Hartree-Fock  determinant  |$o)  in  a canonical  molec- 
ular orbital  basis,  the  negative  of  the  poles  of  G0(E)  are  Koopmans’  values  for  IP’s 
and  EA’s,  which  is  shown  in  the  following.  The  superoperator  algebra  defined  in 
Eq.  3.84  and  Eq.  3.85  is  used.  The  elements  of  the  zeroth  order  inverse  propagator 
matrix  are  given  by 

[Go1],,  = (s\(EI-H0)r) 

= E(s\Ir)  - ( s\H0r ) = £(s|r)  - (s\[r,H0]) 

= £<$0|M]+|$o>  ~ ($o|[r77o,st]+|$o)  + ($0\[H0r,  s']+\<S>0) 

= ^($oksf|$o)  + ^($o|sfr|$o)  - ep($o|rpW|$o) 

p 

- ^2  eP($o|st?’Ptp|$o)  + ep($obVsf|$0)  + ^ eP($o|stpV|$0) 

p p P 

= 5rs(E-er).  (3.110) 
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The  fourth  and  the  fifth  term  do  not  contribute  in  the  line  before  last  in  Eq.  3.110. 
If  r and  s are  occupied  orbitals,  the  first  and  the  third  term  do  not  contribute.  If  r 
and  s are  virtual  orbitals,  the  second  and  the  sixth  term  do  not  contribute. 
Substituting  Eq.  3.106  and  Eq.  3.107  into  Eq.  3.109  gives 

0C(£)  = (Gq1(E)  - E(E))C(E) 

= {(a\(Ei-H0)a)-E(E)]C(E),  (3.111) 


and 

[(a|Jf0a)  + E(£)]C(£)  = EC(E),  (3.112) 

which,  considering  Eq.  3.110,  leads  to 

[e  + Z(E)]C(E)  = EC(E),  (3.113) 


where  e is  a diagonal  matrix  that  contains  the  Hartree-Fock  orbital  energies. 

From  the  matrix  expression  Eq.  3.113,  one  may  abstract  one-electron  equations 
from  the  Fock  operator  / and  the  self-energy  operator  £(E)  [20] 


[/  + E (E)]<pDyson  = E<pDyson.  (3.114) 


The  eigenvectors  are  Dyson  orbitals,  from  which  the  exact  density  matrix  can  be 
constructed  [24] 


D1  Or,*')  = E grf‘DysmW,Er)vD<IKn(x,Er),  (3.115) 

peiP 


where 


dEv{E) 

dE 


E — Ep 


(3.116) 


Since  the  Dyson  orbitals  are  eigenfunctions  of  different  effective  Hamiltonians,  they 
are  not  orthogonal,  and  the  natural  orbitals  (f)p(x)  [66]  and  the  occupation  numbers 
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np  can  be  obtained  by  solving  the  eigenvalue  problem 

J D1(x,x')(f)p(x')dx'  = np(j)v{x).  (3.117) 

If  the  Green’s  function  itself  is  known,  the  ground-state  density  matrix  can  be 
obtained  using  contour  integration  [89],  or  equivalently,  from  sum  rules.  The  contour 
integration  encloses  all  poles  corresponding  to  ionization  potentials 

j G„(E)dE  = 

= <«'o|g,p|4'o)  = D'm.  (3.118) 


Taking  the  first  moment,  leads  to 


1 

2^ n 


EGpq(E)dE  = 


^(Eo-  En)(^\^n)(^n\p\^o) 

n 


+ ^ y^(pr||sf)(^o|gVtts|^0) 

s rst 

= 'y  ] hpsEsq  + - y '(pr\  \st)Dtsrq  = A pq,  (3.119) 

s rst 


where  h is  the  one-electron  part  of  the  Hamiltonian,  given  in  Eq.  2.111.  With 
the  knowledge  of  the  1-matrix  and  the  2-matrix,  the  ground-state  energy  can  be 
calculated 


E0  = Y,h„Dlw  + 1 ^2{rs\\tu)Dlt„  = \tr(hD'  + A).  (3.120) 

pq  rstu 

Coming  back  to  Eq.  3.108,  choosing  the  Hartree-Fock  determinant  to  be  the 
ground-state  wave-function;  a is  then  the  manifold  of  operators  that  produces  the 
configurations  of  the  P space  for  the  IP-EOM  problem  unified  with  the  configurations 
of  the  P space  for  the  EA-EOM  problem.  The  f manifold  corresponds  to  the  unified 
Q spaces  of  the  IP-EOM  and  the  EA-EOM  problem.  Eq.  3.113  can  now  be  compared 
to  Eq.  3.55.  The  self-energy  and  the  correlation  potential  differ  by  the  fact  that  the 
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reference  state  |T0)  in  Eq.  3.113  is  given  by  the  Hartree-Fock  determinant,  whereas 
in  Eq.  3.55  it  is  the  coupled-cluster  wave- function.  In  Eq.  3.113  the  ionization 
potential  and  electron  affinity  components  of  the  electron-propagator  are  treated 
simultaneously,  whereas  Eq.  3.55  is  formulated  for  IP’s,  with  a similar  equation  for 
EA’s,  which  are  decoupled  and  given  in  Eq.  3.51  and  Eq.  3.52. 

3.3  Coupled-Cluster  Green’s  Function  Method 

It  was  shown  by  Meissner  and  Bartlett  [107],  that  the  consequence  of  introduc- 
ing a coupled-cluster  wave-function  as  the  iV-electron  ground  state  in  the  propagator 
approach,  is  the  vanishing  of  the  [G_1(£')]i0  elements,  which  decouples  the  IP  and 
the  EA  sector  for  the  eigenvalues.  The  coupled-cluster  Greens  function  approach 
was  also  discussed  by  Nooijen  and  Snijders  [108,  109]. 

The  reference  function  T0  in  the  propagator  approach  is  chosen  to  be  the 
coupled-cluster  wave-function  T0  = eT<I>0,  with  the  T operator  defined  in  Eq.  3.8 
and  Eq.  3.9.  Using  the  fact  that  strings  of  creation  and  annihilation  operators 
commute,  a new  superoperator  H is  defined 

H = e T HeT  = Ho  + Hcorr,  (3.121) 

where  Hq  is  defined  as  in  Eq.  3.105.  Analogously  to  Eq.  3.25,  the  amplitudes  satisfy 
the  equations 

o = w\mo) 

0 = (3.122) 

With  the  new  defined  superoperator,  Eq.  3.121,  the  inverse  propagator  matrix 
Eq.  3.104  becomes 


G ~\E)  = (a|£/a)  - (a|tfa)  - (a|tff)(f|(EJ  - ^)f)"1(f|^a). 


(3.123) 
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Eq.  3.113  defines  an  effective  eigenvalue  problem 

He//C  (E)  = C (E)E.  (3.124) 

with 

H eff  = (a\Ha)  + ( a\Hi)(t\(EI  - H)i)~\i\Ha).  (3.125) 

For  convenience,  a further  classification  of  the  elements  of  f is  introduced.  The 
two  classes  of  elements  f = {X*,  Y}  are  defined  as 

/ e x*  if  ($0|/  = o 

f e Y if  / |$o)=0.  (3.126) 

The  definition  Eq.  3.126  implies,  that  X*  includes  the  two- hole  one- particle,  three- 
hole  two-particle  and  so  on  operators,  while  Y*  includes  the  two-particle  one-hole, 
three-particle  two-hole  and  so  on  operators. 

Considering  the  [He//]ja  elements  of  the  first  term  in  Eq.  3.125 

(i\Ha)  = ( i\[a,H ]) 

= ($o|[a^,zt]+|$o)  - ($o|[^a,if]+|$o) 

= ($0|a/ftf  + i*aH  — Hai*  — *f^a|$0)  = 0.  (3.127) 

The  first,  the  third,  and  the  fourth  term  in  the  last  line  of  Eq.  3.127  are  zero  because 
of  the  properties  of  the  annihilation  and  creation  operators,  the  second  term  is  zero 
because  of  Eq.  3.122. 

The  first  part  of  the  second  term  of  Eq.  3.125  is  evaluated  for  j^ab  € Y 

(i\Hj]ab)  = (; i\[j]ab,H} ) 

- (<h0|bta^,zt]+|$o)  - <<MtfitaMt]+|$o> 

- ($„|  j'abHi'  + i]j]abH  - Hj'abi'  - i' Hj'ab\%)  = 0.  (3.128) 
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Again,  the  first,  the  third,  and  the  fourth  term  in  the  last  line  of  Eq.  3.129  are 
zero  because  of  the  properties  of  the  annihilation  and  creation  operators,  the  second 
term  is  zero  because  of  Eq.  3.122.  Similarly,  the  three- particle  two-hole  or  higher 
excitation  operators  in  Y give  zero  contributions. 

The  last  part  of  the  second  term  with  a)ij  € X*  is  given  by 

(aUj\Rb)  = (■ aU  j\[b,H ]) 

= ($o\[bHjWa]+\*o)  - ($0\[Hb,jWa]+\$o) 

— (d>0| bHjU^a  + jU^abH  — HbjU^a  — jU^aHb\$o)  = 0.  (3.129) 

The  same  argument  as  above  applies,  as  well  as  for  higher  excitation  operators  in 

xt. 

The  final  part,  needed  to  be  analyzed  to  determine  the  [H e//]ia  elements,  is 

(aUj\Hk^cd)  = (aUj\[k^cd,  H}) 

= ($o\[kfcdH,jHja]+\$o)  ~ (§o\[Hk]cd,  jtita]_H|^»0> 

= ($0|  k^cdHjU^a  + jU^ak^cdH  — Hk^cdjU^a 

—jh^aHtfcd\&0)  = 0.  (3.130) 

Again,  the  only  term  left  in  the  last  line  is  the  second  term,  which  is  zero  because 
of  Eq.  3.122.  Higher  excitation  operators  in  X*  and  Y yield  zero  as  well. 
Summarizing  the  above  results: 

(i\Ra)  = 0 
(X]\Ha)  = 0 
(i\HY)  = 0 
(Xf|^Y)  = 


0. 


(3.131) 
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Eq.  3.131  together  with  the  expansion 


(i\(EI  - H) f)"1  = E-k{£\Ht)k~l 

k= 1 

implies  that 

(i\Hf)(i\(EI  - H)f)~\f\Ha)  = 0, 


(3.132) 


(3.133) 


which  shows  that 

[He//]to  = 0 Vi,  a.  (3.134) 

However,  if  the  same  analysis  for  the  [He//]at  elements  is  done,  it  is  found  that 
the  condition  needed  for  the  [H e//]0i  elements  to  be  zero,  is  that  the  Hermitian 
conjugate  of  Eq.  3.122  is  fulfilled.  Here  only  one  example  is  shown 


(a\Hi)  = (<#,#]) 

= (<h0|[*tf,at]+|<&o>  - <4>o|[Hi,at]+|$0> 

= ($0 VHa]  + aUH  - Hia)  - a] Hi\§ 0).  (3.135) 


The  first,  the  second,  and  the  fourth  term  are  zero  but  for  the  third  term  to  be  zero, 
(<f>o|-H|<f>“)  would  have  to  be  zero.  Since  H is  not  Hermitian,  ( a\Hi ) is  in  general 
nonzero.  The  nonhermiticity  also  results  in 

(a\HX')  ± 0 

(Y| Hi)  ± 0 

(Y|^Xf)  ^ 0.  (3.136) 

It  follows  that 

[He//]ai^0.  (3.137) 

The  [H eff]ij  and  the  [He//]a(,  elements  are  in  general  nonzero.  Hence,  the 
general  structure  of  the  matrix  representation  of  the  effective  Hamiltonian  in  Eq. 
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3.124  is 


(3.138) 


To  obtain  the  eigenvalues  of  the  effective  Hamiltonian,  the  1-hole  sector  can  be 
considered  separately  from  the  1-particle  sector. 

By  introducing  correlation  effects  in  the  ground  state,  the  calculation  of  ioniza- 
tion potentials  and  electron  affinities  decouples.  However,  if  one  is  interested  in  the 
eigenvectors  of  Eq.  3.124,  the  separation  is  not  justified.  For  a complete  decoupling 
of  the  IP  and  EA  sector,  which  also  holds  for  the  eigenvectors,  the  [H eff]ai  and 
the  elements  have  to  be  zero.  In  the  following  the  IP  and  EA  sector  are 

regarded  as  decoupled,  keeping  in  mind,  that  this  is  only  true  for  the  eigenvalues. 

Because  of  Eq.  3.138,  Eq.  3.124  can  be  rewritten  into  an  eigenvalue  equation 
for  the  negative  of  the  IP’s,  and  an  eigenvalue  equation  for  the  negative  of  the  EA’s 


For  the  calculation  of  the  principal  IP’s  the  self-energy  expression  is  simplified 


£/p(£)  - {h\Sh)corr  + {h\HX')(X'\(Ei  - H)X.')-\X^\Hh),  (3.140) 


where  h is  an  array  of  field  operators,  annihilating  occupied  orbitals  in  the  Hartree- 
Fock  determinant.  The  following  shows  that  Eq.  3.140  is  equivalent  to  Eq.  3.54  for 
the  IP  sector. 

Starting  with  the  first  term  in  Eq.  3.140 


[eIP  + XIP(E)]C(E)  = C (E)E 
[eEA  + Y,EA(E)}C(E)  = C {E)E. 


(3.139) 


because  of  Eq.  3.131.  The  manifold  f can  be  replaced  by  the  manifold  X+ 


(k\Hi) 


corr 


-($o| Hitf  + k'Hi  - iHk]  - kUH\$0) 

('hol^  corr  ■ 


corr 


(3.141) 


84 


Further,  assuming  the  truncation  level  of  the  IP-EOM-CCSD  model, 


{k\Ha)ij)  = — (<&0|  HaHjk*  + Hahj  — a\jHk * — k^aUjH  |<f>0) 

= ~{<!>0\k*HaUj\$o)  (3-142) 

(bUk\Hi)  = -($o\HikWb  + tfPbHi-iHtfPb-tffllriH\$0) 


= -{%\kn'bHi\$0). 


(3.143) 


Finally, 

(bUk\HaUj)  = — ($o|  HaHjkH^b  + kH^bHaHj  — aUjHkH^b  — kU*baHjH\$>o) 

= -($0\kU'bHaUj\<S>o).  (3.144) 

Using  Eq.  3.141  through  Eq.  3.144,  the  ij  elements  of  the  self-energy  for  the 
IP  sector  are  given  by 

[EIP(E)]ij  = -(<!>0\ilHj\$0)^r  + ($0\i'HX'\$o)[El 

+($o|X^Xt|$o)]-1($o|X^j|$0>.  (3.145) 

Since  Xf  includes  in  general  the  two-holes  one-particle,  three-holes  two-particles 
and  so  on  operators,  Eq.  3.145  can  be  written  in  matrix  form  with  the  excited 
determinants  of  the  Q space  of  the  IP-EOM  model  as  basis.  Eq.  3.145  can  then  be 
identified  with  Eq.  3.53  for  the  IP  sector,  where  E is  the  negative  of  the  ionization 
potential  e. 

It  was  shown  above,  that  the  [H eff]ai  do  not  vanish  because  of  the  nonhermicity 
of  H.  As  a consequence,  in  particular,  because  of  Eq.  3.136,  the  equations  for 
electron  affinities  are  more  complicated,  since  certain  simplifications  for  the  EA 
sector  do  not  occur  if  H is  used.  However,  the  use  of  H leads  to  the  vanishing  of 
the  [Hg^y]ai  elements,  and  again  a decoupling  of  the  IP  and  EA  sector,  and  similar 
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equations  to  those  for  the  IP’s  can  be  derived  for  the  EA’s.  H 'e^  is  given  by 

K„  = (alA’a)  + (a|Pf)(f|(£/  - *).  (3.146) 

Similar  to  Eq.  3.122,  the  amplitudes  can  be  determined  through  the  following 
set  of  equations 

0 = ($o|£f|$?) 

0 = (^ol^l^6)-.-.  (3.147) 

A similar  analysis,  as  above  for  the  [H eff]ia  elements,  can  be  done  for  the  [H 'efj\ai 
elements. 

As  an  example,  consider 

Hfi’i)  = (<#,#]) 

= - (4>„pV,]+|4'o) 

= ($0 + aW  - WicS  - a]H]i |$0)  = 0.  (3.148) 

The  first,  the  second,  and  the  fourth  term  in  the  last  line  of  Eq.  3.148  are  zero, 
because  of  the  properties  of  the  annihilation  and  creation  operators,  the  third  term 
is  zero,  because  of  Eq.  3.147. 

Proceeding  as  above  will  lead  to 

(altfV)  = 0 
(Y|  h\)  = 0 
(Y|  tfV)  - 


0, 


(3.149) 
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from  which  the  conclusion  can  be  drawn  that  the  [H 'eff]ai  elements  are  zero.  H 'eff 
has  the  following  general  structure 


[K„)n 

[H \ 

0 

[H J 

(3.150) 


That  shows  that  the  decoupling  of  the  IP  and  the  EA  block  is  also  valid  if  H is 
used.  Because  of  Eq.  3.150,  the  self-energy  expression  for  the  EA  sector  can  be 
simplified  as  well 


Vea(E)  = (p|  Afp)corr  + - H]) Y)-1(Y|^tp),  (3.151) 


where  p is  an  array  of  field  operators,  annihilating  virtual  orbitals  in  the  Hartree- 
Fock  determinant. 

Applying  superoperator  algebra  as  above,  leads  to 

[Xea(E)U  = ($o|^tat|$o)corr  + (^o|Y^tat|$o)[^l 

-($o|Y^tYt|$o)]~1(^o|^tYt|<E0).  (3.152) 

Taking  the  Hermitian  adjoint  of  Eq.  3.152,  gives 

[XEA(E)]*ba  = (^o\aHb^\^o)Carr  + {^o\aHY^o)[El 

— ($o|Y^Yt|$o)]_1(^o|Y^’^|$0)-  (3.153) 

Since  Yt  includes  the  two-particles  one-hole,  three-particles  two-holes  and  so  on 
operators,  Eq.  3.153  can  be  written  in  matrix  form,  with  the  excited  determinants 
of  the  Q space  of  the  EA-EOM  model  as  basis.  Assuming  the  correlation  potential 
is  real,  Eq.  3.153  can  then  be  identified  with  Eq.  3.54,  where  E is  the  negative  of 
the  electron  affinity  e. 

It  can  be  concluded,  that  if  the  coupled-cluster  wave-function  is  used  to  de- 
scribe the  N-electron  ground  state,  the  propagator  approach  is  equivalent  to  the 
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partitioned  equation-of-motion  coupled-cluster  approach,  as  far  as  the  eigenvalues 
are  considered. 

3.4  Correlation  Potential  Derived  from  the  Fock  Space 
Coupled-Cluster  Method 

In  the  last  two  sections  a correlation  potential,  Eq.  3.53  and  3.54,  was  derived 
that  yields  the  exact  ionization  potentials  and  electron  affinities  as  the  negative  of  the 
orbital  energies  within  the  correlated  independent  particle  (CIP)  model.  Its  equiv- 
alence to  the  propagator  method  was  shown  if  the  coupled-cluster  wave-function  is 
used  as  the  N-electron  ground  state  and  if  only  the  eigenvalues  are  considered.  How- 
ever, the  derived  correlation  potential  is  energy-dependent  and  consequently,  a dif- 
ferent correlation  potential  corresponds  to  each  orbital,  so  that  the  orbital  equations 
have  to  be  solved  for  each  eigenvalue  separately.  The  eigenvectors  are  nonorthog- 
onal.  In  this  section  a correlation  potential,  using  the  Fock  space  coupled-cluster 
ansatz,  is  derived.  A potential  for  the  IP  sector  and  a potential  for  the  EA  sector 
are  obtained,  both  are  energy- independent  and  universal  for  each  sector.  All  IP’s 
and,  in  an  independent  calculation,  all  EA’s  can  be  obtained  by  a single  diagonal- 
ization  of  an  effective  Hamiltonian.  The  resultant  orbitals  are  biorthogonal,  and  the 
eigenvalues  are  the  negative  of  the  exact  principal  IP’s  and  EA’s. 

In  the  following  the  Fock  space  coupled-cluster  method  is  introduced,  and  it 
is  then  shown,  how  a correlation  potential  can  be  obtained.  Alternatively,  a par- 
titioning technique  is  used  to  obtain  the  Fock  space  coupled-cluster  equations  and 
the  correlation  potential.  In  the  last  part  of  this  section  the  quality  of  the  orbitals 
is  discussed. 

3.4.1  Fock  Space  Coupled-Cluster  Method 

Fock  space  methods,  for  a review  [110],  just  like  equation-of-motion  methods 
[95,  96],  provide  energy  differences  between  electronic  states  directly.  The  Fock 
space  coupled-cluster  method  introduced  by  Lindgren  [92]  and  Haque  and  Mukherjee 
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[111]  is  a multi-reference  method,  in  which  a multiple-cluster  expansion  ansatz  is 
used  to  create  the  different  states  in  Fock  space.  The  Fock  space  is  divided  into 
sectors  (m,  n),  according  to  how  many  electrons  are  added  or  removed  relative  to 
the  reference  function.  We  are  only  interested  in  the  (0, 1)  sector  of  single  ionization 
and  the  (1,0)  sector  of  single  electron  attachment. 

The  model  space  for  the  (N  — 1)  electron  system  contains  a linear  combination 
of  one- hole  determinants,  and  the  model  space  for  the  (N  + 1)  electron  system  is 
a linear  combination  of  one-particle  determinants  relative  to  the  reference  determi- 
nant. Since  we  are  interested  in  all  orbital  energies,  the  number  of  hole  determinants 
is  equal  to  the  number  of  occupied  orbitals  and  the  number  of  particle  determinants 
is  equal  to  the  number  of  virtual  orbitals 


i*n 

occ 

= £cfci«|*0> 

i 

(3.154) 

IO 

vir 

(3.155) 

a 


In  general,  the  Fock  space  CC  method  is  formulated  for  incomplete  model  spaces, 
where  only  certain  hole  and  particle  determinants  are  selected  [110].  The  total 
wave-function  for  the  (N  — 1)  state  and  the  (N  + 1)  state  is  given  by 


l*n  = (3.156) 

IO  = nEA\V°A).  (3.157) 

An  exponential  ansatz  is  used  for  the  wave  operator  Q 

a'p  = {er<«+T<°.»}  (3158) 

aEA  = {erM>+T“’0>}.  (3.159) 

T(°’°)  is  the  normal  closed-shell  coupled-cluster  operator  from  Eq.  3.8.  eT(0,1)  acts  on 


the  model  space  |VF(P)  to  deliver  its  orthogonal  complement  of  net  one- hole  excited 
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determinants,  the  QIP  space  for  the  IP  sector.  er(1,0)  acts  on  the  model  space  I'Pf'4) 
to  deliver  its  orthogonal  complement  of  net  one-particle  excited  determinants,  the 
QEA  space  for  the  EA  sector.  Since  the  (0,  0)  cluster  operator  commutes  with  other 
cluster  operators  in  normal  order,  the  (0,  0)  calculation  can  be  decoupled  from  the 
rest  of  the  Fock  space  calculation 


nIP  = 
nEA  = 


,r<°’°>  .fgT0’1)  j _ eT<0’0)Q/p 

(3.160) 

,r<o,°>  |er(,°)|  _ eTl°’0)£iEA 

(3.161) 

Defining  Hn  — H — Ecc , the  Schroedinger  equation  can  be  written  in  terms  of 
the  model  spaces.  Using  the  model  space  projectors  Pip  and  PEA,  the  Schroedinger 
equation  is  given  by 


HnVLipPip  = nIpuIFPIP  (3.162) 

HNnEAPEA  = QeaujeaPea,  (3.163) 

where  ujip  are  the  principal  IP’s,  and  uEA  are  the  negative  principal  EA’s.  The 
effective  Hamiltonians  Hj,Pj  and  HEAj  can  be  defined  by  projecting  Eq.  3.162  and 
Eq.  3.163  on  the  left  with  PIPQIP  1 and  with  PEA£lEA~1 


PIpnIp-xHNttIPPIP 

= Pip(HnQip)cPip 

= Pipnlp-1nipuippip  = 

Pipuiip  Pip 

TiIP 

Heff 

= pIP(HNnIP)cpIP 

(3.164) 

pEAsiEA~x  hnCieapea 

= PEA(HNttEA)cPEA 

= pEAnEA~1nEAuEApEA 

pEA^jEA  pEA 

ttEA 

Heff 

= pEA(HNnEA)cpEA. 

(3.165) 

To  determine  the  expansion  coefficients  of  U/p  and  £lEA,  Eq.  3.162  and  Eq.  3.163 
are  projected  onto  the  left  by  QIP  and  QEA,  followed  by  insertion  of  Eq.  3.164  and 
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Eq.  3.165,  giving 

QIPHNnIPPIP  - QIPQIPPIPHIepfPIP  = 0 (3.166) 

QEAHNnEAPEA  - QEA(lEAPEAHEAPEA  = 0,  (3.167) 

the  Bloch  equations.  Expanding  the  operators  in  Eq.  3.166  and  Eq.  3.167,  yields 
the  equations  for  the  amplitudes  for  the  (0, 1)  and  for  the  (1,0)  sector.  Both  sectors 
can  be  solved  with  amplitudes  for  the  (0,  0)  sector  alone.  In  Fig.  3-13  and  Fig. 
3-14  the  equations  for  the  effective  Hamiltonians  in  the  IP  and  in  the  EA  sector 
and  the  amplitude  equations  are  given  diagrammatically.  The  interpretation  of  the 
H interaction  vertices  was  given  in  Fig.  3-11.  The  equations  are  written  for  the 
general  case  of  incomplete  model  spaces,  where  only  selected  orbitals  are  active, 
others  remain  inactive.  The  number  of  active  orbitals  is  equal  to  the  number  of  IP’s 
and  EA’s  that  can  be  obtained.  The  closed  arrows  represent  inactive  orbitals,  the 
double  arrows  represent  active  orbitals.  Notice  that  if  all  orbitals  are  chosen  to  be 
active,  the  contribution  from  the  simple  hole  and  particle  amplitudes  is  zero.  Since 
it  is  necessary  for  all  orbitals  to  be  active  to  derive  the  correlation  potential,  further 
consideration  will  involve  only  simple  arrows,  where  it  is  assumed  that  all  orbitals 
are  active. 

If  the  effective  Hamiltonians  are  expressed  in  canonical  Hartree-Fock  orbitals, 
the  leading  contribution  to  the  diagonal  elements  of  — H [pf  and  H EA{  are  the  Hartree- 
Fock  orbital  energies.  That  can  be  understood  by  expanding  the  H interaction  lines 
in  Fig.  3-13  and  Fig.  3-14  according  to  Fig.  3-11,  see  the  discussion  in  Section 
3.1.3.  The  correlation  potential  can  be  found  by  subtracting  the  Hartree-Fock  orbital 
energies  from  the  diagonal  elements  of  the  effective  Hamiltonians 

-HIpfCPP  = (eIP+VIcp)Cpp  = CPPeIpP 
S f//Cpp  = (tEA  + VfA)C  pp  = C pp£EA, 


(3.168) 

(3.169) 
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Figure  3-13.  Equations  for  the  effective  Hamiltonian  and  the  amplitudes  in 
the  IP  sector  for  the  FS-CCSD  method 


where  eIP  and  eEA  contain  again  the  negative  of  the  IP’s  and  the  negative  of  the 
EA’s.  The  correlation  potential  for  the  IP  sector  is  given  in  Fig.  3-15,  and  the 
correlation  potential  for  the  EA  sector  is  given  in  Fig.  3-16. 

The  amplitude  equations  in  Fig.  3-13  and  Fig.  3-14  can  be  solved  iteratively. 
However,  unlike  for  the  EOM  method,  convergence  breaks  down  in  most  cases  due 
to  intruder  states.  Intruder  states  are  eigenstates  that  lie  outside  the  model  space 
but  can  become  quasi-degenerate  with  a reference  space  determinant,  causing  the 
amplitudes  connecting  the  reference  function  with  the  intruder  state  to  become 
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Figure  3-14.  Equations  for  the  effective  Hamiltonian  and  the  amplitudes  in 
the  EA  sector  for  the  FS-CCSD  method 


large.  The  intruder  state  problem  can  be  partially  avoided  for  incomplete  model 
spaces  [112].  Since  it  is  desired  to  obtain  a correlation  potential  for  all  orbitals,  the 
incomplete  model  space  approach,  which  delivers  only  a few  IP’s  or  EA’s  respectively, 
is  of  limited  use  for  us. 

3.4.2  Fock  Space  Coupled-Cluster  Derived  from  a Partitioning  Scheme 

If  the  coupled-cluster  amplitudes,  the  T (0,1\  and  the  T^1,0)  amplitudes  are  given, 
or  somehow  approximated,  Fig.  3-15  and  Fig.  3-16  defines  the  correlation  poten- 
tial. To  avoid  convergence  problems  introduced  by  intruder  states,  the  Fock  space 
coupled-cluster  equations  can  be  cast  into  a different  form. 
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Figure  3-15.  Correlation  potential  in  the  IP  sector  for  the  FS-CCSD  method 


A 


Figure  3-16.  Correlation  potential  in  the  EA  sector  for  the  FS-CCSD  method 


Substitution  of  the  wave  operators  by  SlIP  = eT(0,1)  = 1 + SIP  and  ClEA  = 
)}  = 1 + SEA  in  Eq.  3.164-3.167,  leads  to 


frlP 

Heff 

= PIPH{\  + SIP)PIF 

(3.170) 

0 

= QIPH(  1 + SIP)PIP  - QIPSIPHlPf 

(3.171) 

fjEA 

Heff 

= PEAH(l  + sEA)PEA 

(3.172) 

0 

= QEAH{l  + SEA)PEA-QEASEAHEfA, 

(3.173) 
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where  the  fact  was  used  that  QIPPIP  = QEAPEA  = 0.  To  simplify  the  notation, 
the  subscript  N in  HN  is  omitted.  The  matrix  form  of  Eq.  3.170-Eq.  3.173  is  given 
by 


u IP 
Heff 

— P\IP  1 PfIP  G IP 

— n Pp  npQO 

(3.174) 

0 

— VXIP  1 V*IP  G IP 
~ nQP  + nQQ° 

— S/P(Hpp  + HpgS/P) 

(3.175) 

f tea 

uEA  , ttEAcEA 

— n pp  -|-  ripQO 

(3.176) 

0 

tt E A , xrEAQEA 

~ nQP  + 

qEA/ttEA  . ttEAcEA\ 
a (,rlpp  + rlpQa  J, 

(3.177) 

where  the  definition  of  the  effective  Hamiltonian  is  substituted  into  Eq.  3.175  and 
Eq.  3.177.  If  S/p  is  known,  diagonalization  of  H [pf  gives  all  IP’s,  and  if  SEA  is 
known,  diagonalization  of  HffAf  gives  all  EA’s.  One  could  attempt  to  solve  Eq. 
3.175  and  Eq.  3.177  iteratively  by  taking  as  starting  guesses 

SIP  = (3-178) 

SEA  = (3.179) 

Unfortunately  these  guesses  are  not  sufficient,  and  convergence  cannot  be  achieved 
according  to  my  experience. 

However,  it  is  possible  to  reformulate  the  partitioned  equation-of-motion  eigen- 
value problem  to  find  an  expression  for  S.  It  was  shown  by  diagrammatic  argu- 
ments [113]  and  otherwise  [97]  that  the  IP-EOM-CCSD  and  the  EA-EOM-CCSD 
approaches  are  equivalent  to  the  FS-CCSD  method  for  the  principal  IP’s  or  EA’s, 
solved  in  the  (0, 1)  and  the  (1,0)  sector,  respectively.  In  those  cases  the  exponential 
ansatz  for  the  wave  operator  U reduces  to  a linear  ansatz.  Therefore,  S contains 
only  linear  terms,  and  the  FS-CC  methods  for  IP’s  and  EA’s  can  be  regarded  as  an 
eigenvalue  problem  for  the  transformed  Hamiltonian  H , described  by  Eq.  3.36. 
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The  partitioned  eigenvalue  equation  for  the  (0,1)  sector  is  once  again  used,  but 
slightly  differently, 


u IP  -Clip 

nPp  nPQ 


\ 


V 


-H&  -HIP 


QQ  ) 


( 

\ 


Cpp  C PQ 

CQP  C 


/ 


QQ  ) 


Cpp  C pQ 

CqP  C QQ 


\ ( n \ 

£p  0 


/ 


v 0 £q  / 


, (3.180) 


where  the  equation  was  multiplied  by  minus  one  to  obtain  the  negative  of  the  IP,  e. 
The  two  equations  are  selected  that  depend  upon  £P) 


- HppCpp  - HpqCqp  = CppSp  (3.181) 

— ^ Hq^Cpp  — Hq^Cqp  = C Qp£p.  (3.182) 

Multiplication  of  Eq.  3.181  and  Eq.  3.182  with  Cpp  on  the  right,  leads  to  the 

equivalent  of  Eq.  3.174  and  Eq.  3.175,  where  the  effective  Hamiltonian  yields  the 

negative  of  the  principal  IP’s  upon  diagonalization 


H ^ = -Hpp  — HpgS/p  (3.183) 

0 = -H^p-H^S/p  + S/p(H/ppp  + H^S/f),  (3.184) 

with  the  definitions 

S/p  = CqpCpp  (3.185) 

= CppEpCpp.  (3.186) 

Cqp  and  Cpp-1  can  be  obtained  from  a diagonalization  of  the  transformed 
Hamiltonian  in  the  entire  P and  Q space  for  the  IP  sector.  This  bypasses  any 
explicit  energy  dependence.  Also,  for  a diagonalization  in  the  entire  space,  intruder 
state  problems  do  not  occur,  and  the  effective  Hamiltonian  defined  in  Eq.  3.183  can 
now  be  easily  built. 
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The  partitioned  eigenvalue  equation  for  the  (1,0)  sector  is  given  similarly  to  Eq. 
3.180,  where  the  negative  of  the  EA  is  denoted  by  e 

( 


( XjEA  tt£M  ^ ( 

nPP  nPQ 

mi  mA 


\ QP 


C PP  CpQ 
c QP  C 


QQ  ) 


\ 

( 

£p 

o\ 

/ 

0 

£q  y 

(3.187) 


Notice,  multiplication  by  minus  one  is  not  necessary.  Selecting  the  two  equations 
that  depend  upon  eP  and  multiplying  by  Cpp  on  the  right,  leads  to 


HEA  tt£M  , tt£Mc<.EA 

eft  — Jdpp  -r  Jli  pci o 


leff 


0 = h^  + h^s^ 


PQ'~ 

r EAi 

lQQk 


oEAfrrEA  , u£Ar|£A\ 

o V^ipp  -+-  MpQ»  J, 


(3.188) 

(3.189) 


with  the  definitions 


SEA  = CQPCpp  (3.190) 

Hf//  = CppEpCpp.  (3.191) 

Again  the  Cqp  and  Cpp  1 can  be  obtained  by  diagonalization  of  the  transformed 
Hamiltonian  in  the  entire  P and  Q space  for  the  EA  sector.  Diagonalization  of  H EAf 
yields  the  principal  EA’s. 

From  the  effective  Hamiltonians  defined  in  Eq.  3.183  and  Eq.  3.188,  the  corre- 
lation potentials  can  be  obtained  by  subtracting  Koopman’s  values  as  before,  leading 
to  Eq.  3.168  and  Eq.  3.169,  where  the  effective  Hamiltonians  are  now  built  by  using 
the  eigenvectors  of  the  transformed  Hamiltonian  in  the  P and  Q spaces,  avoiding 
the  difficulty  of  converging  the  iterative  amplitude  equations,  when  all  orbitals  are 
chosen  to  be  active. 

Eq.  3.168  and  Eq.  3.169  define  a correlation  potential  for  all  occupied  or- 
bitals and  a correlation  potential  for  all  unoccupied  orbitals.  Simply  adding  the  two 
contributions  together,  as  schematically  shown  in  Fig.  3-12,  results  in  an  energy- 
independent  correlation  potential  for  the  CIP  model,  which  is  the  same  for  all  orbital 
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energies.  With  Vc  given  in  Fig.  3-12,  the  eigenvalue  problem  is  now  given  by 

fi  e//C  = (e  + Vc)C  = Ce.  (3.192) 

The  matrix  e contains  the  principal  IP’s  and  EA’s  of  the  IP-EOM-CC  and  the  EA- 
EOM-CC  model,  or  equivalently,  the  principal  IP’s  and  EA’s  of  the  FS-CC  method, 
which  are  exact  if  no  truncation  is  introduced. 

As  pointed  out  before,  the  effective  Hamiltonians  in  Eq.  3.51  and  Eq.  3.52, 
which  were  derived  from  the  partitioned  equation-of-motion  approach,  act  in  the 
space  of  N-electron  Hartree-Fock  determinants,  from  which  an  electron  has  been 
removed,  if  the  IP  sector  is  considered,  or  to  which  an  electron  is  added,  if  the  EA 
sector  is  considered.  The  effective  Hamiltonians  in  Eq.  3.168  and  Eq.  3.169  also  act 
in  the  space  of  N-particle  determinants.  However,  as  before,  each  element  of 
and  H is  defined  by 

(3.193) 

(3.194) 

where  the  \<f>i)  are  one-electron  functions,  represented  by  CPP,  solved  for  the  IP 
sector,  and  the  | <j>a)  are  one-electron  functions,  represented  by  CPP,  solved  for  the 
EA  sector.  Again,  Eq.  3.192  can  equivalently  be  regarded  as  an  effective  one- particle 
equation  in  an  orbital  space,  where  the  orbitals  are  solutions  to  the  biorthogonal 
eigenvalue  equation 

HeffHm)  = (/  + vc)\ (f>m)  = £m\<t>m}-  (3.195) 

/ is  the  Fock  operator  and  vc  is  the  correlation  potential,  given  in  matrix  form  in 
Fig.  3-12.  Since  the  effective  Hamiltonian  is  nonhermitian,  diagonalization  yields 
the  right  eigenvectors  | <f>m)  as  well  as  the  left  eigenvectors  (4>m\,  corresponding  to 
the  same  eigenvalue  em.  The  eigenfunctions  can  be  chosen  to  be  biorthonormal 
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(0m|0n)  = Ann-  The  orbital  energies  ((f)m\Heff\(j)m)  = em  are  the  negative  of  the 
principal  IP’s  and  EA’s. 

Since  the  Fock  space  coupled-cluster  method  for  IP’s  and  EA’s  is  equivalent  to 
the  IP-EOM-CC  and  the  EA-EOM-CC  method  respectively,  solving  the  eigenvalue 
equation,  Eq.  3.192,  is  therefore  equivalent  to  solving  Eq.  3.55  for  each  e.  The 
energy  dependence  of  Vc(e)  in  Eq.  3.55  is  removed  by  virtue  of  the  Rayleigh- 
Schrodinger  expansion,  which  is  implicit  in  the  Fock  space  ansatz,  or  by  using  the 
eigenvectors  of  the  full  transformed  Hamiltonian  (in  P and  Q space)  to  build  the 
effective  Hamiltonian.  In  Section  3.3  we  have  also  identified  Vc(e)  in  Eq.  3.55  with 
the  self-energy,  employed  in  propagator  methods  if  one  is  only  concerned  with  the 
eigenvalues.  That  means,  that  the  correlation  potential  in  Eq.  3.192  is  equivalent  to 
the  self-energy,  giving,  if  no  truncation  is  employed,  the  exact  IP’s  and  EA’s.  Since 
the  energy  dependence  is  removed,  all  IP’s  and  EA’s  can  be  obtained  by  a single 
diagonalization,  whereas  the  IP’s  and  EA’s  in  propagator  methods  have  to  each  be 
calculated  separately. 

3.4.3  Orbitals  and  Properties  Obtained  within  the  CIP  Model 

In  this  section  a correlation  potential  for  the  CIP  model  was  derived,  which 
is  energy- independent,  leading  to  a biorthogonal  eigenvalue  equation,  given  by  Eq. 
3.192 

He//C  = (e  + Vc)C  = Ce.  (3.196) 

The  left  and  the  right  eigenvectors  of  the  nonhermitian  matrix  He//  form  a biorthog- 
onal set,  chosen  to  be  orthonormal 

C C = l.  (3.197) 

The  eigenvectors  form  Slater  determinants.  A determinant  |4>)  is  associated  with 
C and  a determinant  ($|  is  associated  with  C.  If  the  effective  Hamiltonian  is  built 
from  a FS-CCSD  ansatz,  the  eigenvalues  are  the  negative  of  the  FS-CCSD  IP’s  and 
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EA’s.  However,  the  question  is  if  the  orbitals  are  improved  as  well.  As  a measure  of 
the  quality  of  the  orbitals  the  calculation  of  the  kinetic  energy  and  dipole  moments 


are  chosen. 

With  the  correlation  potential  given  in  Fig.  3 12,  the  occupied  and  the  unoc- 
cupied blocks  of  the  effective  Hamiltonian  are  decoupled.  However,  Fig.  3-15  and 
Fig.  3-16  show,  that  there  is  an  indirect  contribution  from  the  virtual  space  to 
V£p,  whenever  an  upwards  pointing  line  is  part  of  an  IP  diagram,  and  that  there 
is  an  indirect  contribution  from  the  occupied  space  to  whenever  a downwards 

pointing  line  is  part  of  an  EA  diagram.  Because  of  the  nonhermiticity  of  H eff  its 
eigenvectors  are  biorthogonal,  whereas  Hartree-Fock  orbitals  are  orthogonal,  mean- 
ing that  the  CIP  model  has  two  sets  of  eigenvectors,  whereas  the  Hartree-Fock 
model  has  one  set.  However,  the  eigenvectors  of  the  CIP  model  are  not  expected  to 
be  very  different  from  Hatree-Fock  orbitals,  obtained  by  a unitary  transformation 
of  canonical  Hartree-Fock  orbitals,  since  the  occupied  orbitals  of  the  CIP  model  are 
obtained  by  a rotation  within  the  occupied  space,  and  the  unoccupied  orbitals  of  the 
CIP  model  are  obtained  by  a rotation  within  the  virtual  space  of  the  Hartree-Fock 
orbitals.  Properties  are  therefore  not  expected  to  be  very  different  from  properties 
computed  within  the  Hartree-Fock  model  either.  In  the  following  this  hypothesis 
will  be  tested. 

The  calculation  of  the  kinetic  energy  is  chosen  as  an  example  that  involves  a 
differential  operator.  In  contrast  to  a multiplicative  operator,  its  expectation  value 
depends  on  the  density  matrix  instead  of  the  density.  The  kinetic  energy  operator 
is  given  by 


The  expectation  value  of  the  kinetic  energy  operator  K within  the  CIP  model  can 
be  written  as 


(3.198) 


i=l 


K = ($|jfc|$). 


(3.199) 
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Since  A:  is  a one-electron  operator,  the  kinetic  energy  is 

K = Tr^K),  (3.200) 

where  the  matrices  K and  D1  have  to  be  given  in  a consistent  basis.  For  instance,  if 
K contains  the  kinetic  energy  integrals  in  the  basis  of  the  eigenvectors  of  Eq.  3.196 
(i.e.,  Kmn  = (<j>m\k\(f>n)),  then  D1  is  a diagonal  matrix,  whose  diagonal  elements  are 
the  occupation  numbers  in  the  CIP  model,  two  for  the  occupied  orbitals  and  zero 
otherwise,  where  the  spin  restricted  case  is  assumed.  K and  D1  could  have  also 
both  been  given  in  the  basis  of  atomic  orbitals. 

Table  3-5  lists  the  kinetic  energies  of  several  molecules.  The  geometries  were 
optimized  at  the  MBPT(2)/DZP  level.  The  kinetic  energies  of  the  CIP  model, 
column  II,  are  compared  to  the  Hartree-Fock  kinetic  energies,  column  I,  and  the 
CCSD  kinetic  energies,  column  VI.  To  obtain  the  Hartree-Fock  kinetic  energies,  K 
in  Eq.  3.200  is  given  in  the  basis  of  canonical  Hartree-Fock  orbitals,  in  which  D1  is  a 
diagonal  matrix,  containing  the  occupation  numbers.  The  CCSD  kinetic  energies  are 
calculated  with  the  correlated  density  matrix  in  the  basis  of  molecular  orbitals.  K is 
given  in  the  basis  of  molecular  orbitals  as  well.  Further  comparison  is  made  with  the 
kinetic  energy  obtained  from  the  first  natural  configuration,  column  III.  The  natural 
orbitals  are  calculated  by  diagonalization  of  the  correlated  CCSD  density  matrix. 
The  kinetic  energy  integrals  are  transformed  into  the  basis  of  natural  orbitals,  and 
D1  is  again  a diagonal  matrix,  containing  the  occupation  numbers  two  for  the  first  N 
orbitals  and  zero  otherwise.  Also  included  in  Table  3-5  are  kinetic  energies  computed 
from  a determinant  of  Brueckner  orbitals,  column  IV,  and  kinetic  energies  of  the  CIP 
model,  where  Brueckner  orbitals  are  used  to  formulate  the  effective  Hamiltonian, 
column  V.  For  a more  detailed  discussion  on  Brueckner  orbitals  see  the  following 
chapter. 
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Table  3-5.  Kinetic  energies  in  Hartree,  DZP  basis 


HF 

CIP(HF) 

l.NC 

BR 

CIP(BR) 

CCSD 

HF 

100.03366 

100.03166 

99.98224 

99.89852 

99.87724 

100.11853 

h2o 

75.92108 

75.91744 

75.93392 

75.90440 

75.90061 

76.11443 

nh3 

56.17124 

56.17032 

56.20265 

56.20408 

56.20318 

56.39030 

ch4 

40.14029 

40.14004 

40.16510 

40.17154 

40.17145 

40.33915 

HCN 

92.46336 

92.32461 

92.59806 

92.66512 

92.65981 

92.90719 

Table  3-5  shows  that  the  kinetic  energies  in  the  CIP  model  are  indeed  very- 
similar  to  the  kinetic  energies  obtained  from  the  reference  determinant.  That  is  true 
when  Hartree-Fock  orbitals  are  used  to  build  the  effective  Hamiltonian  (compare 
column  I and  II)  and  also  when  Brueckner  orbitals  are  used  to  build  the  effective 
Hamiltonian  (compare  column  IV  and  V).  A clear  statement,  about  which  determi- 
nant yields  the  best  kinetic  energy,  cannot  be  made  from  the  examples  shown.  The 
Hartree-Fock  determinant  is  slightly  better  for  HF,  for  H20  and  CH4;  the  different 
determinants  yield  similar  kinetic  energies,  and  for  NH3  and  HCN  the  first  natu- 
ral configuration  and  the  Brueckner  determinant  give  better  results  for  the  kinetic 
energy,  compared  to  the  coupled-cluster  kinetic  energy. 

As  a second  example,  the  calculation  of  dipole  moments  was  chosen.  The  clas- 
sical definition  of  the  dipole  moment  is  given  by 

V = '^2qiri,  (3.201) 

i 

where  qi  is  the  charge  with  position  vector  r{.  Quantum  mechanically,  the  electronic 
dipole  operator  is  — J2iLi  ri ■ The  total  dipole  moment  in  the  CIP  model  can  be 
written  as 

N 

ii=-  ZaRa > (3.202) 

i=l  A 

where  the  first  term  is  the  quantum  mechanical  contribution  of  the  electrons  and  the 
second  term  is  the  classical  contribution  of  the  nuclei.  The  three  vector  components 
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Table  3-6.  Dipole  moments  in  Hartree,  DZP  basis 


HF 

CIP(HF) 

l.NC 

BR 

CIP(BR) 

CCSD 

HF 

2.054 

2.054 

1.969 

1.941 

1.941 

1.935 

h2o 

2.102 

2.102 

2.025 

2.010 

2.010 

1.993 

nh3 

1.677 

1.677 

1.649 

1.654 

1.654 

1.631 

HCN 

3.259 

3.259 

3.005 

2.923 

2.923 

2.916 

(x,y,z)  of  the  electronic  contribution  to  the  dipole  moment  are  given  similar  to  Eq. 
3.200,  for  example  for  /ix 

[ix  = TV(DVx),  (3.203) 

where  [ix  contains  the  integrals  of  the  x component  of  the  dipole  moment,  similar 
for  y and  z.  The  classical  contribution  of  the  nuclei  is  added  afterwards. 

The  same  calculations,  as  were  done  for  the  kinetic  energy,  were  performed  to 
obtain  the  components  of  the  dipole  moments,  only  with  different  property  integrals. 
The  results  are  given  in  Table  3-6.  The  example  molecules  have  nonzero  dipole 
moments  only  in  one  direction. 

Table  3-6  shows,  that  the  dipole  moments  in  the  CIP  model  are  the  same 
as  the  ones,  obtained  from  the  reference  determinants.  The  dipole  moments  in 
the  CIP  model,  which  has  the  Hartree-Fock  determinant  as  reference  determinant, 
are  the  same  as  the  Hartree-Fock  dipole  moments  (compare  column  I and  II)  and 
the  dipole  moments  in  the  CIP  model,  which  has  the  Brueckner  determinant  as 
reference  determinant,  are  the  same  as  the  Brueckner  dipole  moments  (compare 
column  IV  and  V).  However,  in  the  case  of  dipole  moments,  it  becomes  apparent 
that  the  dipole  moments  calculated  from  the  first  natural  configuration  and  the 
Brueckner  determinant  are  much  closer  to  the  coupled-cluster  dipole  moments  than 
the  Hartree-Fock  dipole  moments. 
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The  results  for  the  kinetic  energies  and  the  dipole  moments  indicate,  that  the 
orbitals,  or  the  density  matrix  respectively,  in  the  CIP  model  are  only  slightly  modi- 
fied, compared  to  the  reference  orbitals  or  the  reference  density  matrix.  The  kinetic 
energies  calculated  from  the  Hartree-Fock  determinant  seem  to  be  of  about  the  same 
accuracy  as  from  a first  natural  configuration  or  a Brueckner  determinant,  which  is 
not  very  surprising,  since  the  orbitals  of  the  Hartree-Fock  determinant  are  optimized 
with  respect  to  the  total  energy,  which  is  the  sum  of  kinetic  and  potential  energy. 
However,  the  better  choice  over  the  Hartree-Fock  determinant  for  the  calculation  of 
dipole  moments  seems  to  be  a determinant  of  Brueckner  orbitals  or  the  first  natural 
configuration. 


CHAPTER  4 

SECOND-ORDER  CORRELATION  POTENTIAL 

In  the  previous  chapter  possibilities,  of  how  to  obtain  a correlation  potential 
for  the  CIP  model,  which  yields  the  negative  of  the  exact  IP’s  and  EA’s  as  or- 
bital energies,  were  discussed.  The  partitioned  equation-of-motion  approach  led  to 
an  energy- dependent  correlation  potential,  where  an  eigenvalue  equation  has  to  be 
solved  for  each  orbital  with  corresponding  eigenvalue.  If  the  Fock  space  coupled- 
cluster  method  is  employed,  the  energy  dependence  is  removed  by  virtue  of  the 
Rayleigh-Schrodinger  expansion,  which  is  implicit  in  the  Fock  space  ansatz,  or  by 
using  the  eigenvectors  of  the  full  transformed  Hamiltonian  to  build  an  effective 
Hamiltonian. 

However,  in  the  last  section  of  the  previous  chapter,  it  was  pointed  out,  that  even 
though  the  IP’s  and  EA’s  can  be  obtained  at  the  accuracy  of  Fock  space  coupled- 
cluster,  the  orbitals,  and  consequently  the  density  matrix,  in  the  CIP  model  is  not 
significantly  improved,  compared  to  that  of  the  reference  determinant.  That  is, 
because  the  occupied  orbitals  of  the  CIP  model  are  obtained  by  a rotation  of  the 
occupied  orbitals  of  the  reference  function,  and  the  unoccupied  orbitals  of  the  CIP 
model  are  obtained  by  a rotation  of  the  unoccupied  orbitals  of  the  reference  func- 
tion. Properties  computed  with  the  CIP  model  are  therefore  similar  to  properties 
computed  with  the  reference  function. 

The  correlation  potential,  shown  in  Fig.  3 -12,  has  a sector  that  determines 
the  IP’s,  and  an  independent  sector  that  determines  the  EA’s.  This  is  naturally 
the  case,  since  the  FS-CC  methods  in  the  IP  and  in  the  EA  sector,  from  which  the 
correlation  potential  was  derived,  are  two  independent  methods.  The  off-diagonal 
blocks  (the  ia-block  and  the  ai-block)  were  both  set  to  zero.  It  was  shown  for  the 
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IP-EOM  and  the  EA-EOM  methods  (Section  3.3)  that  the  decoupling  is  correct  for 
the  eigenvalues,  since  the  ia-block  is  zero  and  the  ai-block  is  then  arbitrary  (i.e., 
any  choice  of  ai-elements  would  lead  to  the  same  eigenvalues).  Because  the  IP-EOM 
and  the  EA-EOM  methods  are  equivalent  to  the  FS-CC  methods  for  the  IP  and  the 
EA  sector,  the  decoupling  of  the  IP  and  the  EA  sector  in  the  correlation  potential, 
derived  from  the  FS-CC  ansatz  according  to  Fig.  3-12,  is  correct  as  well. 

Because  of  the  independence  of  the  IP  and  the  EA  sector  for  the  FS-CC  method, 
only  the  occupied-occupied  and  the  virtual-virtual  blocks  are  defined  in  the  correla- 
tion potential.  That  was  sufficient  for  the  infinite-order  correlation  potential  of  the 
previous  chapter,  because  the  IP’s  and  the  EA’s  corresponding  to  the  Fock  space 
coupled-cluster  values  could  be  obtained  by  a single  diagonalization.  In  this  chapter 
the  objective  is  to  derive  a finite-order  correlation  potential,  which  depends  on  the 
orbitals,  used  in  an  iterative  procedure,  where  the  orbitals  are  updated  in  each  it- 
eration. For  such  a procedure  it  is  necessary  to  define  the  occupied-virtual  and  the 
virtual-occupied  block  of  the  correlation  potential  as  well,  since  significant  orbital 
changes  are  achieved  by  a rotation  of  the  virtual  into  the  occupied  orbitals.  Apart 
from  the  correctness  of  the  eigenvalues  being  the  negative  of  the  IP’s  and  EA’s,  the 
orbitals  are  desired  to  improve  the  orbitals  of  the  reference  function,  chosen  to  be 
the  Hartree-Fock  determinant.  In  the  last  section  of  the  previous  chapter  results  for 
dipole  moments  were  presented,  which  indicates  that  Brueckner  orbitals  are  better 
than  Hartree-Fock  orbitals  to  calculate  dipole  moments,  which  was  also  found  by 
Hesselmann  and  Jansen  [114].  The  Brueckner  determinant  has  maximum  overlap 
with  the  exact  wave-function,  and  it  can  be  expected  in  general  that  properties,  cal- 
culated from  a Brueckner  determinant,  are  more  accurate  than  properties  calculated 
from  a Hartree-Fock  determinant.  A second  condition  on  the  correlation  potential  is 
therefore  invoked,  which  requires  that  the  resulting  orbitals  are  Brueckner  orbitals. 
The  condition  of  correctness  of  the  IP’s  and  EA’s  and  the  condition  on  the  orbitals 
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are  independent,  since  the  Brillouin-Brueckner  condition  defines  the  partitioning  of 
the  space  in  occupied  and  virtual  orbitals,  whereas  an  arbitrary  rotation  within  the 
occupied  and  within  the  virtual  space  is  of  no  consequence  for  the  partitioning.  In 
this  chapter  the  knowledge  of  the  previous  chapter  is  used  to  define  the  occupied- 
occupied  block  and  the  virtual-virtual  block  of  the  correlation  potential,  and  the 
Brillouin-Brueckner  condition  is  invoked  to  define  the  occupied- virtual  block  of  the 
correlation  potential. 

In  the  first  section  of  this  chapter  the  Brueckner  method  is  summarized,  in 
the  second  section  a second-order  correlation  potential  is  defined  and  in  the  last 
section  it  is  applied  to  calculate  IP’s  and  EA’s  as  well  as  properties  of  some  sample 
molecules. 

4.1  Brueckner  Method 

The  Brueckner  method  originated  in  nuclear  physics  [40]  for  dealing  with  hard- 
core interaction  in  infinite  nuclear  matter.  It  was  introduced  in  quantum  chemistry 
circles  by  Nesbet  [115]  and  by  Lowdin  [41].  Nesbet  required  the  disappearance  of 
single-excited  clusters  from  a trial  wave- function,  which  he  called  the  Brueckner 
condition.  Lowdin  assumed  the  best-overlap  condition  to  the  exact  wave-function 
in  his  ’’exact  self-consistent  theory”,  resulting  in  the  vanishing  of  single-excitation 
contributions,  which  he  called  the  Brillouin-Brueckner  theorem. 

The  Brillouin-Brueckner  condition  was  used  as  a starting  point  in  the  derivation 
of  orbital  equations,  based  on  different  types  of  correlated  wave-functions.  Larsson 
[42,  43]  used  pair  functions  to  obtain  Hartree-Fock  like  equations,  but  with  an  inho- 
mogeneity term  added.  Kvasnicka  [116]  used  perturbation  theory  to  obtain  Brueck- 
ner orbitals.  Chiles  and  Dykstra  [44],  Purvis  and  Bartlett  [45]  and  Adamowicz  and 
Bartlett  [46],  and  Handy  et  al.  [47]  obtained  Brueckner  orbitals  in  coupled-cluster 
theory.  Stolarczyk  and  Monkhorst  [48,  117,  118]  applied  the  maximum  overlap 
condition  to  a coupled-cluster  reference  wave-function  and  formulated  an  effective 
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Hamiltonian,  which  has  the  form  of  the  Fock  operator  plus  a correlation  potential. 
Scuseria  [49,  50]  implemented  the  Brueckner  coupled-cluster  method,  employing 
the  effective  Brueckner  Hamiltonian  in  a self  consistent  scheme.  Further  aspects  of 
Brueckner  theory  were  recently  presented  by  Lindgren  and  Salomonson  [51,  119], 
including  its  approximate  relationship  to  DFT. 

The  Brillouin- Brueckner  condition  is  only  sufficient  to  determine  the  partition 
of  the  full  space  into  occupied  and  virtual  parts.  For  a complete  specification  of  the 
corresponding  eigenvalue  problem,  the  one-particle  operator  also  has  to  be  defined 
in  the  occupied-occupied  block  and  in  the  virtual- virtual  block  by  applying  an  ad- 
ditional condition.  Scuseria  [49,  50]  required  that  the  total  coupled-cluster  energy  is 
reproduced  as  a sum  of  orbital  energies,  whereas  Kvasnicka  [116]  and  Lindgren  [51] 
assumed  that  the  orbital  energies  correspond  to  IP’s  and  EA’s. 

The  Brueckner  single  determinant  |<f>B)  is  defined  by  requiring  that  the  exact 
wave- function  |T)  has  maximum  overlap  with  |$s) 

(T  |$B)  = max.  (4.1) 

If  a reference  determinant  |<F0)  is  chosen,  then  a determinant,  different  from  |<L0),  is 
obtained  by  rotating  the  virtual  orbitals  into  the  occupied  orbitals,  represented  by 
the  operator  c“{ aU },  if  only  linear  terms  are  considered.  A general  single  determi- 
nant is  then  given  by 

|*)  = (l  + ef{o*i})|«0>.  (4.2) 

For  the  overlap  to  be  maximal,  it  has  to  be  stationary  with  respect  to  the  rotation 

r\ 

— <^|(1  + c“{af*})  |<l?o)  = 0 Vi,  a,  (4.3) 

which  leads  to 


<#!«?>  = 0 Vi,  a. 


(4.4) 
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If  the  orbital  rotation  is  allowed  to  contribute  beyond  linear  terms,  then,  ac- 
cording to  the  Thouless  parameterization  [120],  any  two  nonorthogonal  Slater  de- 
terminants, within  a given  basis,  are  related  by  the  transformation 


with  Ti  given  in  Eq.  3.9.  For  the  overlap  to  be  maximal,  all  7\  amplitudes  have  to 


based  Brueckner  method. 

In  the  Brueckner  CC  method  the  Tx  amplitudes  can  be  rotated  away  iteratively  [45, 
46].  Solving  first  the  coupled-cluster  equations,  the  amplitudes  are  used  to  trans- 
form the  orbitals,  which  give  new  integrals,  with  which  a new  coupled-cluster  cal- 
culation is  performed,  yielding  new  7\  amplitudes.  The  resulting  7\  amplitudes  are 
used  again  to  transform  the  orbitals.  The  procedure  is  repeated  until  the  Ti  ampli- 
tudes reach  a cut  off.  The  orbitals  of  the  last  transformation  are  then  the  Brueckner 
orbitals. 

Alternatively,  a Brueckner  effective  Hamiltonian  can  be  defined,  as  for  instance 
by  Scuseria  [49,  50],  who  constructed  an  effective  Hamiltonian  by  requiring  that  the 
coupled-cluster  energy  can  be  reproduced  as  a sum  of  expectation  values  of  single- 
particle operators.  If  the  truncation  level  of  the  coupled-cluster  wave-function  is 
chosen  to  be  T = T\  + T2,  and  assuming  the  Brueckner  condition  Tx  = 0,  the  total 
energy  is  the  sum  of  the  reference  energy  and  the  correlation  energy,  Eq.  3.23, 


|!<)  = e7'  |<I>o), 


(4.5) 


vanish,  meaning  that  all  effects  of  eTl  must  vanish  as  well,  which  is  the  difference 
between  the  Brueckner  coupled-cluster  method  and  the  configuration  interaction 


E — Eq  + AEccsd  — 


]C(*w*} + \ 5>iiu>  + \j2(ii\\ab)t' 


ab 
ij  ‘ 


(4.6) 


ijab 


The  total  energy  can  be  rewritten  as 


(4.7) 
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v — _x 


Figure  4-1.  Occupied-occupied  block  of  the  brueckner  Hamiltonian 

A 

A x 


Figure  4-2.  Virtual-virtual  block  of  the  brueckner  Hamiltonian 


which  defines  the  Brueckner  effective  Hamiltonian  for  the  occupied-occupied  block, 

Fa  = fij  + \ (4.8) 

kab 

where  fij  is  the  Fock  matrix  element.  Particle-hole  symmetry  arguments  lead  then 
to  the  definition  of  the  effective  Hamiltonian  in  the  virtual-virtual  block 

= (4.9) 

ijc 

The  occupied-virtual  block  is  chosen  to  be  the  T\  equation,  Fig.  3-5,  where  terms 
involving  the  7\  amplitudes  are  zero 

ft . = /*. + E + \ - \ <4-10) 

jb  jbc  jkb 

Convergence  of  the  Brueckner  coupled-cluster  equations  therefore  enforces  = 
Fia  = 0.  The  different  blocks  of  the  Brueckner  Hamiltonian  are  diagrammatically 
given  in  Fig.  4-1,  in  Fig.  4-2,  and  in  Fig.  4-3.  The  equation  for  the  T2  amplitudes 
were  given  diagrammatically  in  Fig.  3-6,  Fig.  3-7,  and  Fig.  3-8,  where  now  all 
diagrams  containing  a Ti  amplitude  are  zero. 

A calculation  can  be  performed  by  starting  with  a set  of  T2  amplitudes,  from 
which  the  Brueckner  Hamiltonian  is  constructed.  It  is  then  diagonalized  to  obtain 
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Figure  4-3.  Occupied-virtual  block  of  the  brueckner  Hamiltonian 


Figure  4-4.  Effective  Hamiltonian  in  the  IP  sector  for  the  FS-CCSD  method 
with  T\  = 0 

a new  set  of  orbitals,  which  define  new  integrals,  which  are  used  to  calculate  new  T2 
amplitudes.  The  procedure  is  repeated  until  self-consistency. 

However,  instead  of  requiring  that  the  couple-cluster  energy  can  be  repro- 
duced,one  can  choose  to  employ  an  effective  Hamiltonian,  which  yields  the  IP’s  and 
EA’s  as  the  negative  of  the  eigenvalues.  A correlation  potential  for  the  occupied- 
occupied  block,  Fig.  3-15,  and  for  the  virtual-virtual  block,  Fig.  3-16,  were  defined 
in  the  previous  chapter,  which  does  just  that,  if  the  Hartree-Fock  terms  are  added. 
If  the  Brueckner  condition  is  superimposed  (Tj  = 0),  the  diagrams  for  the  occupied- 
occupied  and  the  virtual- virtual  block  of  the  effective  Hamiltonian  are  given  in  Fig. 
4-4  and  Fig.  4-5,  which  contain  Fock  space  amplitudes.  Comparing  those  to  Fig. 
4-1  and  to  Fig.  4-2  shows  that  the  difference  lies  in  the  last  two  diagrams  in  Fig. 
4-4  and  Fig.  4-5.  If  the  coupled-cluster  amplitudes  and  the  Fock  space  amplitudes 
are  given,  the  effective  Hamiltonian,  defined  by  Fig.  4-3,  Fig.  4-4,  and  Fig.  4-5, 
yields  Brueckner  orbitals  and  exact  FS-CC  IP’s  and  EA’s  as  eigenvalues. 


Ill 
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EA 

eff 
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\ 

Figure  4-5.  Effective  Hamiltonian  in  the  EA  sector  for  the  FS-CCSD  method 
with  Ti  = 0 


4.2  Second-Order  Brueckner  Hamiltonian  Yielding  Second-Order  IP’s 

and  EA’s 

In  this  section  the  Brueckner  Hamiltonian  yielding  FS-CC  IP’s  and  EA’s  is 
expressed  up  to  second  order  in  perturbation,  which  yields  second-order  Brueckner 
orbitals  as  eigenfunctions  and  IP’s  and  EA’s,  which  are  correct  through  second 
order  as  eigenvalues.  Lindgren  and  Salomonson  [51,  119]  derived  a second-order 
Brueckner  potential,  yielding  correct  second-order  IP’s  and  EA’s,  by  minimizing  the 
energy  expectation  value  with  respect  to  the  partitioning  of  the  Hamiltonian. 

In  order  to  derive  a second-order  Brueckner  potential,  the  Hamiltonian  has  to 
be  partitioned  into  an  unperturbed  part  and  a perturbation.  The  Hamiltonian  in 
second  quantization  was  given  in  Eq.  2.111 

H = ^2hpqP^Q+\^2{PQ\\rs)pfqisr.  (4.11) 

pq  pqrs 

A symmetric  one-electron  operator 

U = ^2uPqp^q  (4.12) 

pq 

can  be  added  and  subtracted  to  give 

H = (Hx  + U)  - U + H2  = F - U + H2, 


(4.13) 
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where  Hi  is  the  one-electron  part  and  H2  is  the  two-electron  part  of  the  Hamiltonian 
in  Eq.  4.11.  F is  given  by 


F Hi  + U 'y  S^PQ  T upq)P^ Q ~ 'y  ] fpqP^Qi 

pq  pq 

(4.14) 

where  / = h + u.  u is  arbitrary  but  chosen  to  be 

U y ^(^i  ^i)i 

i 

(4.15) 

so  that 

UPQ  = ^(piWqi)- 
i 

(4.16) 

The  coulomb  and  exchange  operators  are  defined  by 

Mmi)  = (*(2)i— i*(2))0(i) 

r 12 

(4.17) 

Kt(  1)0(1)  = <*(2)|  — |<A(2))*(1), 

n 2 

(4.18) 

where  it  is  integrated  over  the  coordinates  of  electron  2 only.  This  definition  does  not 
restrict  the  orbitals  to  be  Hartree-Fock  orbitals.  However,  if  Hartree-Fock  orbitals 
are  used,  this  choice  makes  the  operator  F equal  to  the  Fock  operator,  in  which  case 
(canonical)  fpq  = 0 for  p ^ q.  In  general  F is  nondiagonal,  with  fpq  ^ 0. 

In  order  to  define  the  unperturbed  Hamiltonian  H0,  F is  split  into  diagonal  and 
off-diagonal  parts 

F = Fd  + F°.  (4.19) 

The  diagonal  and  off-diagonal  parts  are  given  by 

pi  = £/>*!>  (4-20) 

P 

F°  = E f>'i- 

P^Q 


(4.21) 
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With  these  definitions  the  Hamiltonian  takes  the  form 

H = E fppP]p  + E%  “ um)P^  + \ E<wl I rs)p]q]sr  = H0  + V,  (4.22) 

p pq  pqrs 

where 

Ho  = Fd  = Y,idrrp'v  (4.23) 

V 

V = F0-U  + H2  = '^(f°g-upq)pfq  + ^J2(pqHrs)pfqfsr.  (4.24) 

pq  pqrs 

Normal-ordered  operators  can  be  introduced  by  applying  Wick’s  theorem,  Sec- 
tion 3.1.1,  which  leads  for  the  zero-order  part  to 

Ho  = E tfpptP  = E tfp&P}  + E £l 

P p i 

(H„)n  = if„-B(°)  = ^/4{ptp},  (4.25) 

P 

where  e*.  The  perturbation  part,  Eq.  4.24,  gives 

v = Ew  ~ upq){p](i)  - E uh  + \ E<wiira>w«*} 

pq  i pqrs 

+ eE^ii^hpM + ^ E<viiv> 

pq  i ij 

v = FZ-Un  + Wn  + Un  + ($  o\V\$o) 

V = F%  + Wn+  ($0|V|$o),  (4.26) 

where 


f°n 

= E fnip'fi 

pq 

(4.27) 

UN 

= Euw(pt«}  = EE<p<ii9<){pt9} 

pq  pq  i 

(4.28) 

WN 

= ^^2{PQ\\rs){ptqtsr} 

pqrs 

(4.29) 

<$o|V|$0) 

= -E«-  + 5E<w>- 

(4.30) 

i ij 
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The  normal-ordered  perturbation  operator  is  then 

VN  = V - <$0|^|^o) 

- F°N  + wN 

= ^fpq{p*q}  + \^2{pq\\rs){pWsr}.  (4.31) 

p^q  pqrs 

The  particular  partitioning  of  the  Hamiltonian,  defined  in  Eq.  4.25  and  Eq. 
4.31,  is  referred  to  as  partitioning  I.  Diagrammatically  the  operators  FN  and  WN 
were  given  in  Fig.  3-1  and  Fig.  3-2.  In  order  to  express  the  effective  Hamiltonian, 
defined  in  Fig.  4-3,  Fig.  4-4,  and  Fig.  4-5,  up  to  second-order  in  perturbation,  all 
diagrams  have  to  be  collected  that  contain  up  to  two  vertices.  Expressing  the  Tj 
equation  up  to  second  order,  ensures  that  the  sum  t[x+2}  = T1(1)-)-T1(2)  = 0.  The  indi- 
vidual terms  are  not  necessarily  zero,  and  terms  containing  Tj^  have  to  be  included, 
even  though  terms  that  contain  Tj  disappear  in  the  infinite-order  expressions. 

The  virtual-occupied  block  of  the  second-order  effective  Hamiltonian  for  parti- 
tioning I is  given  in  Fig.  4 6.  The  horizontal  lines  symbolize  the  denominators.  The 
algebraic  interpretation  of  the  diagrams  can  be  found  in  Table  4-1.  For  closed-shell 
systems,  the  spin-restricted  form  can  be  used,  which  is  obtained  by  spin  integration. 
Since  the  aa  effective  Hamiltonian  is  equivalent  to  the  P/3  effective  Hamiltonian, 
only  one  of  them  needs  to  be  constructed.  The  spin-restricted  expressions  can  be 
found  in  Table  4-1  as  well,  where  capital  letters  refer  to  spatial-orbital  indices.  Note 
that  the  occupied- virtual  block  is  defined  as  the  adjoint  of  the  virtual-occupied  block. 

The  occupied-occupied  and  virtual-virtual  block  of  the  effective  Hamiltonian 
are  derived  by  expanding  the  amplitudes  in  Fig.  3-15  and  Fig.  3-16  in  orders 
of  perturbation.  The  terms  up  to  second  order  are  identified  as  the  Fock  matrix 
elements  and  all  diagrams  that  contain  two  interaction  vertices.  The  resulting  di- 
agrams are  given  in  Fig.  4 7 and  Fig.  4 8.  Table  4-2  and  Table  4-3  contain  the 
algebraic  expressions  of  the  diagrams.  With  the  exception  of  the  diagrams  I and  III 
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V VI  VII 

Figure  4-6.  Occupied- virtual  block  of  the  second-order  effective  Hamilto- 
nian, partitioning  I 

in  Fig.  4-7  and  Fig.  4-8,  all  other  diagrams  have  a nonhermitian  contribution.  In 
order  to  obtain  a hermitian  effective  Hamiltonian,  each  term  is  multiplied  by  |,  and 
the  corresponding  term,  where  i and  j and  a and  b respectively  are  exchanged,  is 
added.  Table  4-2  and  Table  4-3  also  show  the  spin-integrated  expressions. 

The  second-order  Brueckner  Hamiltonian  according  to  partitioning  I,  defined  in 
Fig.  4-6,  Fig.  4-7,  and  Fig.  4-8,  is  the  same  as  the  effective  Hamiltonian  Lindgren 
and  Salomonson  [51,  119]  derived  by  minimizing  the  energy  expectation  value  with 
respect  to  the  partitioning  of  the  Hamiltonian. 

However,  if  the  Hamiltonian  is  partitioned  according  to  Eq.  4.25  and  Eq.  4.31 
(partitioning  I),  then  ( Ho)n  is  not  invariant  with  respect  to  a transformation  of  the 
orbitals  within  the  occupied  or  the  virtual  space;  a necessity  for  any  consistent  theory 
of  electron  correlation.  If  the  Hamiltonian  is  partitioned,  so  that  the  occupied- 
occupied  and  the  virtual-virtual  off-diagonal  Fock  matrix  elements  are  part  of  the 
unperturbed  Hamiltonian  (H0)N,  and  not  V/v,  as  in  partitioning  I,  ( H0)N  is  invariant 
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Table  4-1.  Algebraic  interpretation  of  the  diagrams  in  Fig.  4-6 


spin  orbital  expressions 

Ik 


spatial  orbital  expressions 


I 

II 

III 

IV 

V 

VI 


fii~  fa 


fabfbi 

b^a  (fii—fbb){fii~faa) 

Ei'^j 


fa  fa 


E* 


jV*  {fii  faa){fjj~faa) 

(ja\\bi)fbj 

jb  {fjj  fbb){fii  faa) 

(ab\\ij)fjb 

jb  {fii~\~  fjj  ~ faa~  fbb){fii~  faa) 


E 


-IE 


(jk\\ib){ab\\jk) 


2 jkb  ( fii~\~fkk  faa  fbb){fii  faa) 


vn  | E?- 


(qj||fec)<&c||tj) 


2 (fii+fjj  fbb  fee)  {fii  ~ faa) 


Ja  I 


Sii-Saa 


E 


/ab/bj 


B±A  (///-/bb)(///-/aa) 


-E 


SjiSaj 


E 


(///-/aa)(/jj-/aa) 

(2(JA\BI)-(JA\IB))fBJ 


E, 


(/jj— /bb)(/b- /aa) 


'JB  {fn+fjj-fAA-fBB){flI~fAA) 

~2  E JBTB 

(2(JiC|BJ)  — (JA'|/B))(AB|A'.7) 


(2(jrsri/B)-(JBr|BJ))(/iB|jjt') 

{flI+fKK-fAA-fBB){flI-fAA) 


+ 


{flI+fKK-fAA-fBB){flI-fAA) 


IE, 


(2{2U|BC)-(AJ|CB))(BC1/J) 
if II  A"  fjj  f BB  fcc){f 1 1 f Aa) 


+ 


(2pU|CB)-(A.7[BC))(BC|J/) 

{flI+fjJ-fBB~fcc){flI-fAA) 


to  an  orbital  transformation  within  those  subspaces.  Partitioning  II  is  given  by 

<//„)„  = f*  = E + E /£,{j>'9}+  E £,{*}  <4-32) 

p p,g£occ  p,q£vir 

vn  = ^2fpq{piQ}  + ^^2(pq\\rs){piqhr}.  (4.33) 

p,q  pqrs 

The  prime  in  the  sum  means,  that  if  p is  an  index  of  an  occupied  orbital,  then  q is 
an  index  of  a virtual  orbital,  and  if  p is  an  index  of  a virtual  orbital,  then  q is  an 
index  of  an  occupied  orbital. 

Since  the  off-diagonal  fy  and  fab  elements  are  included  to  all  orders  in  the 
unperturbed  Hamiltonian  in  partitioning  II,  the  counting  of  orders  becomes  slightly 
different.  Diagram  II  in  Fig.  4-6,  diagram  III  in  Fig.  4-7  and  Fig.  4-8  are  first- 
order  diagrams  for  partitioning  II,  instead  of  second-order  diagrams  for  partitioning 
I.  Partitioning  II  also  implies,  that  a first-order  in  Vjv  solution  for  the  amplitudes, 
which  includes  the  off-diagonal  /y  and  fab  elements  to  infinite-order,  has  to  be 


117 


V VI  VII 

Figure  4-7.  Occupied-occupied  block  of  the  second-order  effective  Hamilto- 
nian, partitioning  I 

found,  before  the  second-order  terms  can  be  built.  Fig.  4-9  shows  the  iterative 
equations  for  the  first-order  single  excitation  and  double  excitation  amplitudes.  The 
algebraic  interpretation  can  be  found  in  Table  4-4,  together  with  the  spin-integrated 
expressions.  For  the  amplitudes,  the  permutation  of  i and  j in  diagram  V and 
a and  b in  diagram  VI  has  to  be  taken  into  account.  Table  4-4  therefore  shows  two 
terms  for  diagram  V and  two  terms  for  diagram  VI.  The  spin  integration  of  the  T[ 1 > 
amplitudes  is  similar  to  the  spin  integration  of  the  effective  Hamiltonian,  only  the  aa 
part  needs  to  be  considered.  The  spin  integration  of  the  doubly  excited  amplitudes 
leads  to  two  different  types  of  amplitudes:  aaaa,  denoted  by  T^\  and  aa(3/3, 
denoted  by  [T^]',  which  are  iterated  separately.  For  simplicity,  the  superscript  (1) 
is  omitted  in  the  algebraic  expressions. 

Fig.  4-10  shows  the  virtual-occupied  block  of  the  second-order  effective  Hamil- 
tonian, employing  partitioning  II  in  terms  of  the  first-order  amplitudes.  The  occupied- 
virtual  block  is  again  its  adjoint.  The  interpretation  of  the  diagrams  in  spin-orbital 
and  spin-integrated  form  is  given  in  Table  4 5. 
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Figure  4-8.  Virtual-virtual  block  of  the  second-order  effective  Hamiltonian, 
partitioning  I 


In  principle,  the  first-order  Fock  space  amplitudes  should  be  iterated  to  include 
all  zeroth-order  terms,  as  has  been  done  for  the  coupled-cluster  amplitudes.  How- 
ever, since  the  does  not  change  significantly  during  the  iteration,  the  first-order 
Fock  space  amplitudes  are  simply  given  as  in  Fig.  4-11,  with  the  algebraic  inter- 
pretation in  Table  4-6.  The  spin  integration  leads,  as  for  the  amplitudes,  to 
two  different  types  of  Fock  space  amplitudes:  denotes  the  aaaa  spin  case  and 

[Si1*]'  denotes  the  aa(3[3  spin  case. 

Fig.  4-12  shows  the  occupied-occupied  block  of  the  second-order  Hamiltonian 
according  to  partitioning  II,  Fig.  4-13  shows  the  virtual- virtual  block.  The  algebraic 
interpretations  are  given  in  Table  4-7  and  Table  4-8  in  the  spin-orbital  and  the 
spin-integrated  from.  In  order  to  obtain  an  hermitian  effective  Hamiltonian,  each 
diagram,  except  diagram  I,  in  Fig.  4-12  and  Fig.  4-13  is  multiplied  by  a factor  of 
one  half,  and  the  diagram,  in  which  i and  j and  a and  b respectively  is  exchanged, 
is  added. 
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Table  4-2.  Algebraic  interpretation  of  the  diagrams  in  Fig.  4-7 


spin  orbital  expressions 


spatial  orbital  expressions 


I 

II 

III 

IV 

V 

VI 

VII 


fi 


Ea 


y 

fig  ft 


aj_ 


fjj  fa 


^2k^i,j  f ikfkj 

'ka  fkk—fa 


E(ik\\ja)fak 
k 


(ia\\jk)fka 


^—'ka  fkk+fjj—faa—fii 

1 X^  (fc/||ja>(ia||fcO 

2 


2 2—lkla  fkk+fll  — faa—fii 


1 X~^  <fci||6a><6ffl||fcj) 

2 Z-/kab  fkk+fjj-faa-fbb 


flJ 

Efl  A f AJ 
A fjj  f AA 

~ f IK  f KJ 

V (2(IK\JA)-(IK\AJ))fAK 


;KA 


f KK  f AA 


y-  (2{IA\JK)-(IA\KJ))fKA 

2-^KA  f kkA~ f j J f AA  f II 


(2(KL\JA)-{KL\AJ))(IA\KL) 

fKK+fLL-fAA-fll 


"2  E KLA 

(: 2{KL\AJ)-{KL\JA)){IA\LK ) 


+ 


£E 


KAB 


f KkA-  f LL  f AA  f 11 


(2(KI\BA)-(KI\AB))(BA\KJ) 

f KK  “T  f J J f AA  f BB 


+ 


(2(KI\AB)-{KI\BA))(BA\JK) 

f K K f J J ~ f AA~  f B B 


Depending  on  the  choice  of  the  partitioning  of  the  Hamiltonian,  two  second- 
order  effective  Hamiltonians  were  derived.  Its  eigenfunctions  are  second-order  Brueck- 
ner  orbitals  and  its  eigenvalues  are  the  negative  of  the  ionization  potentials  and 
electron  affinities,  correct  through  second-order 

H^C<2>  = C<V2)  (4.34) 

The  effective  Hamiltonian  depends  on  the  orbitals,  which  leads  to  an  iterative 
scheme.  Starting  with  Hartree-Fock  orbitals,  the  effective  Hamiltonian  is  built  and 
diagonalized.  The  new  orbitals  are  used  to  build  new  integrals  and  a new  Fock 
matrix,  with  which  a new  effective  Hamiltonian  is  built.  Its  eigenfunctions  are 
again  used  to  transform  the  integrals  and  to  build  a Fock  matrix.  The  procedure  is 
repeated  until  self-consistency. 
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Table  4-3.  Algebraic  interpretation  of  the  diagrams  in  Fig.  4-8 


I 

II 

III 

IV 

V 

VI 


spin  orbital  expressions 


fab 

Efaifib 

facfcb 

E{ai\\bc)fci 

(ac\\bi)fcj 
fii~\~fbb~faa~fcc 

1 (ai\\cd)  (cd\\bi) 

2 2—jicd  fu+fbb—fcc—fdd 


£ 


spatial  orbital  expressions 
Jab 

_ V"  S AlflB 

Sii-Saa 

£c^A,B  f AcfcB 

v (2(AI\BC)  — (AI\CB))fci 
Z-'IC  f ii -fee 

y*  (2(AC\BI)  — (AC\IB))Jci 

fll+fBB-fAA-feC 

1 ^ [ (2{AI\CD)  — (AI\DC))(CD\BI) 

2 A-'ICD  " fll+fBB-fce-jDD 

_j_  (2(AI\DC)—(AI\CD))(CD\IB) 
f 1 1 f b b — fee  f dd 


VII 


1 04|co)  (cb\\ji) 

2 A^ijc  fii+fjj—fcc  — fbb 


1 (2{JI\CA)—(JI\AC))(CB\JI) 

2 2-iIJC  fnA f j j — fee — f b b 

+ (2(J  I\AC)  — (J  I\C  A)){C  B\1  J)  1 
f ii  A fjj  fee — f b b 


4.3  Numerical  Results 

In  the  previous  section  it  was  shown,  how  to  obtain  second-order  Brueckner 
Hamiltonians,  which  yield  second-order  IP’s  and  EA’s,  according  to  two  different 
partitioning  schemes.  For  practical  reasons,  the  second-order  CIP  model  is  im- 
plemented slightly  differently.  Since  the  Brueckner  condition  defines  the  orbitals, 
and  diagonalization  of  the  occupied-occupied  and  the  virtual- virtual  block  is  merely 
a unitary  transformation  of  the  orbitals,  the  procedure  is  implemented  in  a two- 
step  scheme.  First,  the  second-order  Brueckner  orbitals  are  determined  iteratively, 
second,  the  occupied-occupied  and  the  virtual-virtual  blocks  of  the  effective  Hamil- 
tonians are  built  with  the  previously  obtained  orbitals.  Diagonalization  yields  then 
the  second-order  IP’s  and  EA’s. 

For  partitioning  I the  orbital  determination  consist  of  the  following  steps: 

1.  transformation  of  the  one-electron  integrals 

2.  transformation  of  the  two-electron  integrals 
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Figure  4-9.  Iterative  equations  of  the  first-order  single  and  double  excitation 
amplitudes,  partitioning  II 


3.  construction  of  the  Fock  matrix 

4.  construction  of  the  second-order  7\(2)  amplitudes,  as  given  in  Fig.  4-6  and  Table 
4 1 

5.  rotation  of  the  orbitals 

6.  check  if  the  orbitals  are  converged;  if  yes,  stop;  if  no,  go  back  to  1. 

For  partitioning  II  the  orbital  determination  is  slightly  different,  since  the  first- 
order  amplitudes  have  to  be  determined  iteratively  as  well: 

1.  transformation  of  the  one-electron  integrals 

2.  transformation  of  the  two-electron  integrals 

3.  construction  of  the  Fock  matrix 

4.  iteration  of  the  first-order  amplitudes,  according  to  Fig.  4-9  and  Table  4-4 
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Table  4-4.  Algebraic  interpretation  of  the  diagrams  in  Fig.  4-9 


I 

II 

III 

IV 

V 


spin  orbital  expressions 

fii~  faa 

Efajfji  J.CL 
fii~faa  3 

Efabfbi  j-b 
fii—faa  3 

(ab\\ij) 

fjj  faa  fbb 


E 


M. 


fiiA  fjj— faa— fbb  ^ 


ab 


Y'  . hi fab 

fii+Jjj—faa  — fbb  fk 


spatial  orbital  expressions 
Sai 

ill  f AA 

Ef aj i j i +A 

I^I  f 1 1 ~ f AA  J 

ESabSbi  iB 
B*A  Su-fAAlJ 

1).  (AB\IJ)-(AB\JI) 

2 ’ fll  A-  fjj  f AA  f BB 

fr{1)i/-  (ab\ij) 

1 2 1 ' f II A f J J f A A~  f BB 

V itu. +ab 

2 ’ fuAfjJ-fAA-fBB  IK 

lb1’]'  = t2(i) 

Lk± +ab 

2 flJ+fjJ-fAA-fBB  JK 

Pa1’]'  = 0 


VI  E 


fbc 


fii  A fjj  — faa  fbb  0 

Y''  fac f 

L^ic^a  fiiAfjj-faa-fbb  t? 


rpW.TT' 

J2  -Z^ 

* T2(1):-E 


f BC 


4-AC 


C^B  fu+fjJ—fAA—fBB  II 

[t2(1)]'  = t2(1) 


f AC 


-t 


BC 


C¥^A  fuAfjJ-fAA-fBB  II 

[T2(1)]'  = 0 


5.  iteration  of  the  first-order  X 2^  amplitudes  for  the  aaaa  spin  case,  according  to 
Fig.  4-9  and  Table  4-4 

6.  iteration  of  the  first-order  [T2(1)]'  amplitudes  for  the  aaf3/3  spin  case,  according  to 
Fig.  4-9  and  Table  4-4 

7.  construction  of  the  second-order  amplitudes,  as  given  in  Fig.  4-10  and  Table 
4-5 

8.  rotation  of  the  orbitals 

9.  check  if  the  orbitals  are  converged;  if  yes,  stop;  if  no,  go  back  to  1. 
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I II  III 


IV  V 

Figure  4-10.  Virtual-occupied  block  of  the  second-order  effective  Hamilto- 
nian, partitioning  II 


The  orbital  rotation  step  is  carried  out  in  such  a way,  that  the  new  occupied 
orbitals  i'  are  orthogonal  to  the  new  virtual  orbitals  a' 


10  = N)  + X^ala> 

a 

la'>  = (4-35) 

i 

However,  the  occupied  orbitals  and  the  virtual  orbitals  are  not  orthogonal  among 
themselves.  To  obtain  orthogonal  orbitals  m",  the  orbitals  n'  are  transformed,  ac- 
cording to 

\m")  =^TXnm\n'),  (4.36) 

n 

where  X = S~2,  and  S is  the  overlap  matrix  S = (C/)t  C'.  It  can  be  easily  verified, 
that  this  transformation  indeed  orthogonalizes  the  orbitals 


(C")f  C" 

s-^c7)1  c's-5 


(crs-iVcs-* 

S-5SS-5  = 1. 


1 

2 


(4.37) 
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Table  4-5.  Algebraic  interpretation  of  the  diagrams  in  Fig.  4-10 


spin  orbital  expressions 


spatial  orbital  expressions 


II 


I 


(ja||&t)  .b 
'-—‘jb  fu-faaLj 


III 


IV 


1 v MMl+cfe 

2 l^jbc  fa-fan  Lij 


V 


1 (jk\\ib)  ,ab 

2 Z^jkb  fa  faa  jk 


As  a measurement  of  the  quality  of  the  orbitals,  the  calculation  of  the  dipole 


molecules  are  given.  For  comparison,  the  dipole  moments  obtained  from  the  Hartree- 
Fock,  the  MBPT(2),  and  the  CCSD  density  matrices  are  included,  together  with  the 
results,  obtained  from  the  second-order  Brueckner  orbitals,  when  partitioning  I and 
partitioning  II  is  used,  as  well  as  the  results  after  the  first  step  of  the  iteration, 
which  is  the  same  for  partitioning  I and  II,  indicated  by  nonit.  The  DZP  basis  set 
was  used  for  all  calculations.  The  geometry  of  the  molecules  was  optimized  at  the 
MBPT(2)/DZP  level  of  theory. 

Fig.  4-14  shows  the  differences  of  the  CCSD  dipole  moments  to  the  dipole 
moments,  obtained  from  the  other  methods,  included  in  Table  4-9.  The  dipole 
moments,  obtained  from  the  Brueckner  determinant,  are  closest  to  the  CCSD  re- 
sults. The  second-order  Brueckner  approximations  are  a significant  improvement 
over  Hartree-Fock,  and  close,  but  somewhat  worse,  than  the  MBPT(2)  results.  This 
is  not  surprising,  since  in  the  MBPT(2)  method  the  dipole  moments  are  calcu- 
lated from  a second-order  correlated  density  matrix,  whereas  in  the  second-order 


moment  is  chosen,  Section  3.4.3.  In  Table  4-9  the  dipole  moments  of  some  sample 
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II 

Figure  4-11.  first-order  Fock  space  amplitudes,  partitioning  II 


Table  4-6.  Algebraic  interpretation  of  the  diagrams  in  Fig.  4-11 


spin  orbital  expressions 

spatial  orbital  expressions 

I 

(ia\\jk) 

cHT  (IA\JK)-(IA\KJ) 

fkk~faa~fii 

2 ' f JJ+fKK-f AA-fll 

rc(i)i/.  (ia\kj) 

L 2 1 ‘ f JJ+fKK-f AA~f II 

II 

(ac\\bi) 

c(l).  (AC\BI)-{AC\IB) 

fbb~  faa~  fee 

2 ' fll+fBB-fAA-fcC 

roU)y.  (AC\Bi) 

2 * ’ fll+fBB-fAA-fcC 

Brueckner  methods,  the  1-particle  functions  are  approximated  up  to  second  order 
in  perturbation,  and  the  dipole  moments  are  calculated  from  an  idempotent  density 
matrix.  The  calculation  of  the  dipole  moment  of  CO  is  a particular  difficult  case, 
and  the  errors  are  big  compared  to  the  other  test  cases;  but  all  approximations 
correctly  change  the  sign  of  the  dipole  moment,  which  is  predicted  wrong  by  the 
Hartree-Fock  method.  In  order  to  view  the  differences  between  the  second-order 
Brueckner  approximations  better,  in  Fig.  4-15  the  dipole  moment  differences  are 
given,  where  CO  is  taken  out,  and  the  Hartree-Fock  errors  for  HCN  and  H2CO  are 
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Figure  4-12.  Occupied-occupied  block  of  the  second-order  effective  Hamilto- 
nian, partitioning  II 


cut  off  at  0.15  Deb.  The  graph  shows,  that  the  second-order  Brueckner  orbitals, 
obtained  by  calculating  them  from  Hartree-Fock  orbitals,  according  to  Fig.  4-6  or 
Fig.  4-10,  are  close  to  the  iterative  solutions.  Partitioning  II  yields  usually  better 
results  than  partitioning  I. 

After  the  orbitals  are  determined,  the  occupied-occupied  and  the  virtual- virtual 
blocks  of  the  second-order  Brueckner  Hamiltonians  are  built  and  diagonalized  to 
obtain  the  negative  of  the  second-order  IP’s  and  EA’s  as  eigenvalues.  Table  4-10, 
Table  4-11,  Table  4-12,  Table  4-13,  Table  4-14,  Table  4-15,  and  Table  4-16  show 
the  eigenvalues  of  the  second-order  Brueckner  Hamiltonians,  according  to  partition- 
ing I and  II  for  some  sample  molecules.  The  IP’s  and  EA’s  in  the  tables  are  given  in 
Hartree  to  avoid  numbers  with  large  magnitude.  Also  included  in  the  tables  are  the 
eigenvalues  when  Hartree-Fock  orbitals  are  used  to  build  the  second-order  Hamilto- 
nian (nonit.),  which  is  identical  for  the  two  partitioning  schemes.  The  eigenvalues  are 
compared  to  Hartree-Fock  values,  as  well  as  IP-EOM-CCSD  and  EA-EOM-CCSD 
results.  The  infinite-order  EA-EOM-CCSD  results  are  chosen  as  reference  values  in 
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Table  4-7.  Algebraic  interpretation  of  the  diagrams  in  Fig.  4-12 


spin  orbital  expressions 

spatial  orbital  expressions 

I 

fij 

fu 

II 

Ea 

Ea  fAjtj 

III 

EkamjaK 

T,ka  (1(1K\JA)  - (IK\AJ))t& 

IV 

Elea  fkas]ak 

EKA  f KA  ( Sjk  + [sJkY ) 

V 

2 EkAB  [(KI\BA)  {pKJ  + [^k'j]') 

+{KI\AB)  (t%  + [(“]')] 

VI 

-lEkla(kl\\ja)Skl 

~\  Y.KLA  \(KL\JA)  (SffL  + [Sj^L]') 

+ (KL\AJ)  (s™  + [sft]')] 

the  L 2 space  for  the  second-order  approximations,  even  though  we  don’t  claim  that 
those  unbound  states  exist  (some  might  in  the  continuum).  The  DZP  basis  set  was 
used,  and  the  geometries  were  optimized  at  the  MBPT(2)/DZP  level  of  theory. 

In  Fig.  4-16,  Fig.  4-18,  Fig.  4-20  and  Fig.  4-22  the  differences  in  eV  between 
the  negative  of  the  IP-EOM-CCSD  results  for  the  IP’s  and  the  eigenvalues  of  the 
different  independent  particle  models,  included  in  the  tables,  for  HF,  H20,  NH3,  and 
CH4  (Table  4-10,  Table  4—11,  Table  4-12,  Table  4-13)  are  shown.  The  core  IP’s 
are  not  given,  since  the  energy  differences  are  relatively  large,  because  of  their  large 
absolute  values.  All  energy  differences  in  the  graphs  are  presented  in  eV  (IP’s  are 
commonly  reported  in  eV).  Fig.  4 17,  Fig.  4-19,  Fig.  4-21,  and  Fig.  4-23  show  the 
differences  in  eV  between  EA-EOM-CCSD  results  for  the  EA’s  and  the  eigenvalues  of 
the  independent  particle  models.  Each  cluster  of  columns  belongs  to  a particular  IP 
or  EA,  for  which  the  EOM-CCSD  result  in  eV  is  given  next  to  the  cluster.  The  first 
observation  to  be  made,  is  that  the  second-order  Hamiltonians  generally  overcorrect 
Koopmans’  values  for  IP’s  and  EA’s.  While  the  eigenvalues  of  the  occupied  orbitals 
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Figure  4-13.  Virtual-virtual  block  of  the  second-order  effective  Hamiltonian, 
partitioning  II 


are  too  low  and  the  eigenvalues  of  the  virtual  orbitals  are  too  high  for  the  Hartree- 
Fock  model,  the  second-order  Hamiltonians  predict  eigenvalues  for  the  occupied 
orbitals  that  are  too  high,  and  for  the  virtual  orbitals  that  are  too  low.  The  second- 
order  approximations  for  IP’s  and  EA’s  are  often  very  close.  However,  for  the  EA’s 
of  H20  at  5.3  eV  and  21.5  eV  partitioning  I shows  particular  large  deviations  from 
the  EA-EOM-CCSD  results.  A large  deviation  from  the  EA-EOM-CCSD  result 
occurs  for  the  noniterative  scheme  for  H20  at  21.5  eV.  The  values  for  IP’s  and 
EA’s,  obtained  from  partitioning  scheme  II,  are  of  relatively  stable  quality  and  in 
most  cases  an  improvement  over  the  Hartree-Fock  values.  In  general,  the  EA’s  are 
better  described  by  the  second-order  approximation  than  the  IP’s,  compared  to 
the  Hartree-Fock  values,  which  is  the  result  of  the  error  cancellation  occurring  for 
Koopmans’  IP’s. 

However,  the  second-order  IP’s  and  EA’s  for  CO,  H2CO,  and  HCN,  Table 
4-14,  Table  4-15,  and  Table  4-16,  are  of  rather  poor  quality.  It  becomes  difficult 
to  correlate  the  second-order  values  with  the  EOM-CCSD  results,  even  the  wrong 
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Table  4-8.  Algebraic  interpretation  of  the  diagrams  in  Fig.  4-13 


spin  orbital  expressions 

spatial  orbital  expressions 

I 

fab 

Jab 

II 

- Ei  M 

- E / flBtf 

III 

EiCH \bc)tCi 

T.,c(2{AI\BC)  - (AI\CB))tf 

IV 

E icficSbi 

£,cM4?  + K?]') 

V 

-5Eyc0*l|a*>$ 
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VI 
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IE, CD  [{AI\CD)  (.gf  + [«gfl') 

+ (A1\DC)  (sf°  + [ogf]')] 

signs  are  predicted.  The  difference  between  those  molecules  and  the  ones  discussed 
before,  is  that  they  all  contain  multiple  bonds,  which  means  the  probability  of 
intruder  state  problems  increases.  Intruder  state  problems  are  common  in  Fock 
space  coupled-cluster  and  were  discussed  briefly  in  Section  3.4.1.  For  instance,  with 
a standard  Fock  space  coupled-cluster  code  only  3 roots  of  H20  can  be  converged 
in  a DZP  basis  set.  The  second-order  approximations  seem  to  be  less  sensitive  to 
intruder  state  problems,  but  they  still  are  apparent. 

Having  a closer  look  at  the  IP’s  of  CO,  Table  4-14  column  3,  one  can  see,  that 
the  two  eigenvalues  with  value  -1.30494  H for  partitioning  II  do  not  seem  to  fit. 
Assuming  that  those  are  intruder  states,  leaving  them  out,  and  changing  the  order 
of  the  two  eigenvalues  at  -0.58524  H and  the  eigenvalue  at  -0.57059  H,  a correlation 
between  eigenvalues  of  the  different  methods  and  the  EOM-CCSD  results  can  be 
found.  Fig.  4-24  shows  for  the  six  lowest  IP’s  the  difference  in  eV  between  the  IP- 
EOM-CCSD  values  and  the  other  approximate  methods.  The  experimental  order 
of  the  lowest  states  confirms  the  IP-EOM-CCSD  results  (5cr  14.01  eV,  l7r  16.91  eV, 
4cr  19.72  eV)  [121].  The  usual  trend  is  found,  that  the  second-order  approximations 
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overcorrect  the  Hartree-Fock  values.  Focusing  on  the  partitioning  scheme  II,  the 
overcorrection  leads  for  the  IP’s  at  19.4  eV  and  13.9  eV  to  a larger  difference  to 
IP-EOM-CCSD  than  Hartree-Fock. 

The  same  can  be  done  for  the  IP’s  of  H2CO,  Table  4-15  column  3.  Leaving 
the  eigenvalue  for  partitioning  II  at  -1.59545  H out,  leads  to  the  energy  differences 
in  eV,  shown  in  Fig.  4-25.  The  lowest  IP-EOM-CCSD  results  are  confirmed  by 
experiment  (b2  10.9  eV,  bi  14.5  eV,  bi  16.1  eV,  b2  17.0  eV,  bi  21.4  eV)  [121].  Again, 
the  Hartree-Fock  values  are  generally  overcorrected,  with  larger  differences  to  IP- 
EOM-CCSD  for  the  second-order  approximations  at  10.4  eV,  14.2  eV,  and  15.7  eV 
than  Hartree-Fock. 

For  HCN,  Table  4-16,  it  is  very  difficult  to  see  which  eigenvalues  belong  to 
which  IP  and  it  was  not  attempted  to  correlate  the  different  states. 

Fig.  4-24  and  Fig.  4-25  show,  that  when  intruder  states  are  taken  into  ac- 
count, it  is  possible  to  make  the  connection  between  the  EOM-CCSD  values  and  the 
second-order  approximations.  However,  the  intruder  states  also  influence  the  other 
eigenvalues  through  the  diagonalization  of  the  Hamiltonians  [122].  Included  in  all 
the  tables  are  the  percentages  of  singles.  If  the  percentage  is  high,  the  corresponding 
IP  root  has  mainly  one-hole  character.  The  smaller  the  percentage  gets,  the  more 
two- hole  one-particle  character  the  root  has  (i.e.  an  intruder  state  gains  influence), 
which  lowers  the  quality  of  the  root.  A high  percentage  of  singles  for  EA’s  means, 
that  the  character  of  the  root  is  mainly  one-particle  character,  a low  percentage 
indicates  an  increasing  influence  of  a two-particle  one-hole  state. 

In  Table  4-10,  Table  4-12,  and  Table  4-13  the  percentage  of  singles  is  high  for  all 
roots  of  HF,  NH3,  and  CH4.  In  relation,  the  second-order  IP’s  and  EA’s  consistently 
improve  Koopmans’  values.  For  H20  (Table  4-11)  the  percentage  of  singles  for  the 
second  highest  IP’s  is  only  60.4%.  Consequently,  inconsistencies  in  the  second-order 
eigenvalues,  especially  for  the  EA’s,  occur.  For  CO  the  third  highest  IP  has  only  a 
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Table  4-9.  Dipole  moments  in  Deb 


HF 

part.  I 

part.  II 

nonit. 

Brueckner 

MBPT(2) 

CCSD 

HF 

2.054 

1.948 

1.937 

1.947 

1.941 

1.945 

1.935 

h2o 

2.102 

2.023 

2.013 

2.026 

2.010 

2.012 

1.993 

nh3 

1.677 

1.665 

1.661 

1.664 

1.654 

1.655 

1.631 

HCN 

3.259 

2.937 

2.901 

2.939 

2.923 

2.894 

2.916 

CO 

0.369 

-0.439 

-0.376 

-0.357 

-0.188 

-0.310 

-0.064 

H2CO 

2.934 

2.310 

2.285 

2.353 

2.366 

2.326 

2.387 

value  of  53.2%  singles,  for  H2CO  that  value  is  43.8%,  and  for  HCN  43.8%  (Table 
4-14,  Table  4-15,  Table  4-16)  and  significant  errors  occur  in  the  IP  and  the  EA 
spectra. 
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Figure  4-14.  Dipole  moment  differences  from  CCSD  dipole  moments  in  Deb 
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Figure  4-15.  Dipole  moment  differences  from  CCSD  dipole  moments  in  Deb, 
the  Hartree-Fock  values  for  HCN  and  H2C0  are  cut  off 
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Table  4-10.  Negative  IP’s  and  EA’s  in  Hartree  for  HF/DZP 


HF 

part.  I 

part.  II 

nonit. 

IP/EA-EOM-CCSD 

%S 

-26.28539 

-25.05174 

-25.08629 

-25.08724 

-25.65810 

85.7 

-1.59133 

-1.39889 

-1.41129 

-1.40721 

-1.43576 

82.9 

-0.75860 

-0.66569 

-0.67669 

-0.67275 

-0.71583 

95.7 

-0.64098 

-0.50330 

-0.51584 

-0.51080 

-0.56560 

94.9 

-0.64098 

-0.50330 

-0.51584 

-0.51080 

-0.56560 

94.9 

0.22143 

0.20110 

0.19894 

0.20041 

0.20349 

98.7 

0.98659 

0.91893 

0.91275 

0.91538 

0.92105 

94.5 

1.05931 

1.01886 

1.01039 

1.01504 

1.01912 

96.7 

1.05931 

1.01886 

1.01039 

1.01504 

1.01912 

96.7 

1.14869 

1.10746 

1.10342 

1.10536 

1.10492 

95.2 

1.50856 

1.47320 

1.47204 

1.47277 
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Figure  4-16.  IP  differences  from  IP-EOM-CCSD  in  eV  for  HF/DZP 
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Figure  4-17.  EA  differences  from  EA-EOM-CCSD  in  eV  for  HF/DZP 


Table  4-11.  Negative  IP’s  and  EA’s  in  Hartree  for  H20/DZP 


HF 

part.  I 

part.  II 

nonit. 

IP/EA-EOM-CCSD 

%S 

-20.55751 

-19.43178 

-19.44382 

-19.45110 

-19.96103 

84.2 

-1.33196 

-1.15820 

-1.16357 

-1.16282 

-1.17460 

60.4 

-0.71173 

-0.65495 

-0.65948 

-0.65786 

-0.68814 

95.8 

-0.56036 

-0.46204 

-0.46725 

-0.46518 

-0.51003 

94.7 

-0.49922 

-0.38784 

-0.39324 

-0.39024 

-0.43616 

94.4 

0.22131 

-0.04088 

0.16231 

0.18866 

0.19516 

97.9 

0.31264 

0.28087 

0.27742 

0.27987 

0.28496 

97.4 
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0.72253 

0.77832 

0.73692 

0.78943 
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0.86854 
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0.86318 

0.85997 

0.86132 

0.86787 

90.4 

1.02785 

0.96681 

0.96572 

0.96528 

0.97146 

92.0 

136 


> 

<D 


2.00 

1.00 

0.00 


-1.00 


-2.00 


-3.00 


-4.00 


-5.00 


2a,  32.0 


0 


I 


I 


lb,  18.7 


i 


i 


3a,  13.9 


i 


i 


lb2  11.9 


□ Koopmans’  values 
H Partitioning  I 
EJ  Partitioning  II 
0 noniterative 


Figure  4-18.  IP  differences  from  IP-EOM-CCSD  in  eV  for  H20/DZP 


Table  4-12.  Negative  IP’s  and  EA’s  in  Hartree  for  NH3/DZP 


HF 

part.  I 

part.  II 

nonit. 

IP/EA-EOM-CCSD 

%S 

-15.54241 

-14.64534 

-14.63836 

-14.64866 

-15.00636 

83.3 

-1.13735 

-1.00089 

-0.99994 

-1.00092 

-1.01459 

81.3 

-0.62069 
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-0.57474 
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95.4 
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-0.57474 
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-0.59681 

95.4 
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-0.35657 
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Figure  4-19.  EA  differences  from  EA-EOM-CCSD  in  eV  for  H20/DZP,  value 
for  partitioning  I at  5.3  eV  is  cut  off 


Table  4-13.  Negative  IP’s  and  EA’s  in  Hartree  for  CH4/DZP 


HF 

part.  I 

part.  II 

nonit. 

IP/EA-EOM-CCSD 

%S 

-11.21049 

-10.57459 

-10.56399 

-10.57498 

-10.76674 

82.8 

-0.94181 

-0.84581 

-0.84338 

-0.84531 
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Figure  4-20.  IP  differences  from  IP-EOM-CCSD  in  eV  for  NH3/DZP 
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Figure  4-21.  EA  differences  from  EA-EOM-CCSD  in  eV  for  NH3/DZP 
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Table  4-14.  Negative  IP’s  and  EA’s  in  Hartree  for  CO/DZP 


HF 

part.  I 

part.  II 

nonit. 

IP /EA-EOM-CCSD 

%S 

-20.67026 

-19.33478 

-19.40594 

-19.44733 

-20.07916 

84.1 

-11.38191 

-10.91112 

-10.85325 

-10.88541 

-10.98422 

83.2 

-1.50623 

-1.22295 

-1.30494 
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-1.31468 
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Figure  4-22.  IP  differences  from  IP-EOM-CCSD  in  eV  for  CH4/DZP 
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Figure  4-23.  EA  differences  from  EA-EOM-CCSD  in  eV  for  CH4/DZP 
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Figure  4-24.  IP  differences  from  IP-EOM-CCSD  in  eV  for  CO/DZP 
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Table  4-T5.  Negative  IP’s  and  EA’s  in  Hartree  for  H2CO/DZP 


HF 

part.  I 

part.  II 

nonit. 

IP/EA-EOM-CCSD 

%S 

-20.57847 

-19.24650 

-19.29043 

-19.32706 

-19.97270 

83.3 

-11.35315 

-10.78168 

-10.74734 

-10.76450 

-10.91576 

82.3 

-1.39698 

-1.20105 

-1.59545 

-1.10494 

-1.17432 

43.8 

-0.86952 

-0.77683 

-1.14507 

-0.77595 

-0.79056 

88.5 

-0.68874 

-0.60810 

-0.76830 

-0.60098 

-0.62874 

91.3 
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-0.49467 
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Figure  4-25.  IP  differences  from  IP-EOM-CCSD  in  eV  for  H2CO/DZP 
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Table  4-16.  Negative  IP’s  and  EA’s  in  Hartree  for  HCN/DZP 


HF 

part.  I 

part.  II 

nonit. 

IP/EA-EOM-CCSD 

%S 

-15.61928 

-14.57737 

-84.58809 

-14.60710 

-15.07147 

82.1 

-11.31547 

-10.69301 

-14.57316 

-10.68916 

-10.88826 

82.8 

-1.22330 

-8.91337 

-10.69232 

-5.25436 

-1.05901 

47.4 

-0.81278 

-5.26210 

-10.33742 

-0.70557 

-0.74736 

90.7 

-0.58113 

-0.71036 

-0.69820 

-0.48857 

-0.50258 

91.1 

-0.48710 

-0.49220 

-0.48794 

-0.48857 

-0.49657 

95.7 

-0.48710 

-0.49220 

-0.48794 

0.08215 

-0.49657 

95.7 

0.15607 

0.08370 

-0.30914 

0.08215 

0.13216 

95.5 

0.15607 

0.08370 

0.08208 

0.13958 

0.13216 

95.5 

0.22862 

0.16644 

0.08208 

0.35395 

0.19493 

96.6 

0.41458 

0.32914 

0.27179 

0.38467 

0.38976 

96.7 

0.47164 

0.37629 

0.33715 

0.38467 

0.42234 

93.6 

0.47164 

0.37629 

0.33715 

0.42320 

0.42234 

93.6 

CHAPTER  5 
CONCLUSION 

In  this  dissertation  the  correlated  independent  particle  (CIP)  model  is  discussed. 
The  CIP  model  can  be  regarded  as  a generalized  Hartree-Fock  model,  which  includes 
correlation,  or  as  an  approximation  to  density  matrix  functional  theory. 

A correlated  independent  particle  (CIP)  model  is  formulated,  which  has  the 
general  form  (/  + vc)(j)p  = Ep<fip,  where  / is  the  Fock  operator  and  vc  is  the  corre- 
lation potential  to  be  determined,  while  the  4>p  are  the  one-particle  solutions.  The 
first  question  to  be  answered  is,  what  conditions  should  be  imposed  leading  to  the 
construction  of  the  correlation  potential.  That  condition  is  chosen  to  be  the  ex- 
actness of  the  ionization  potentials  and  electron  affinities,  which  correspond  to  the 
eigenvalues  of  the  model.  This  choice  is  inspired  and  supported  by  the  extended 
Koopman’s  theorem. 

In  Chapter  2,  following  the  introduction,  the  basic  principles  and  ideas,  lead- 
ing to  the  formulation  of  the  CIP  model,  are  discussed.  Theorems  and  proofs  are 
summarized,  which  are  essential  to  density  matrix  functional  theory  (DMFT).  In 
DMFT  the  total  energy  is  written  as  a functional  of  the  density  matrix,  and  the 
energy  is  minimized  with  respect  to  the  density  matrix.  Even  though  DMFT  is 
similar  in  spirit  to  density  functional  theory,  where  local  potentials  are  employed, 
the  potentials  used  in  DMFT  are  general,  giving  the  possibility  for  nonlocal  poten- 
tials. While  the  kinetic  energy  and  the  exchange  term  are  known  in  DMFT,  the 
correlation  energy  functional  is  unknown. 

The  extended  Koopmans’  theorem  is  summarized  in  Chapter  2 as  well.  While 
Koopmans’  theorem  provides  a simple  one-electron  model  for  ionization  and  electron 
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attachment,  the  extended  Koopmans’  theorem  is  a generalization,  which  uses  a cor- 
related wave-function  to  describe  the  unionized  or  unattached  system.  The  extended 
Koopmans’  theorem  employs  a nonlocal  and  nonhermitian  effective  one-particle  op- 
erator, whose  eigenvalues  are  the  variationally  optimized  energy  differences  between 
the  ionic  and  parent  species,  i.e.  the  ionization  potentials  and  the  negative  of  the 
electron  affinities. 

In  Chapter  3 the  condition  of  exactness  of  the  ionization  potentials  and  elec- 
tron affinities  as  eigenvalues  is  exploited  to  obtain  correlation  potentials.  If  no 
truncation  is  introduced,  the  equation-of-motion  coupled-cluster  method  for  ion- 
ization potentials  and  electron  affinities  (IP-EOM-CC/EA-EOM-CC)  yields  exact 
ionization  potentials  and  electron  affinities.  An  energy-dependent  correlation  poten- 
tial is  obtained  from  a partitioned  equation-of-motion  approach,  which  reproduces 
IP-EOM-CC  and  EA-EOM-CC  results  as  orbital  energies.  As  a consequence  of  the 
energy  dependence  of  the  correlation  potential,  a different  effective  Hamiltonian  cor- 
responds to  each  orbital  and  each  eigenvalue.  An  eigenvalue  equation  has  therefore 
to  be  solved  for  each  eigenvalue  separately.  The  resulting  orbitals  are  nonorthogo- 
nal,  and  the  density  matrix  is  not  idempotent,  even  though  the  general  form  of  the 
CIP  model  as  shown  above  is  maintained. 

The  energy  dependence  is  regarded  as  a disadvantage,  but  can  be  removed  by 
using  the  Fock  space  coupled-cluster  method  to  extract  a correlation  potential.  The 
Fock  space  coupled-cluster  method  for  the  one-hole  and  one-particle  sector  and  the 
IP-EOM-CC/EA-EOM-CC  method  are  equivalent,  and  the  eigenvalues  of  the  effec- 
tive Hamiltonians  are  the  same,  if  the  same  truncation  level  is  used.  The  correlation 
potential  is  constructed  from  the  Fock  space  coupled-cluster  amplitudes.  Since  those 
are  difficult  to  converge  for  all  ionization  potentials  and  electron  affinities,  an  alter- 
native route  to  build  the  correlation  potential  is  shown,  which  uses  the  eigenvectors 
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of  the  transformed  Hamiltonian  in  the  full  configuration  space.  The  resultant  cor- 
relation potential  is  energy-independent  and  universal  for  all  orbitals.  In  a one-step 
procedure  the  exact  FS-CC  results  for  ionization  potentials  and  electron  affinities 
are  obtained  within  the  CIP  model,  using  the  correlation  potential  derived  from  the 
Fock  space  coupled-cluster  method. 

Also,  in  Chapter  3 the  electron-propagator  method  is  discussed.  In  contrast 
to  the  EOM-CC  and  Fock  space  methods,  the  one-hole  and  the  one-particle  sectors 
are  in  general  not  decoupled,  but  described  simultaneously.  From  the  electron- 
propagator  method  a formally  exact  one-electron  theory  can  be  derived,  where  an 
energy-dependent  effective  Hamiltonian  is  employed,  which  is  the  sum  of  the  Fock 
operator  and  the  energy-dependent  self-energy.  It  is  shown,  that  if  the  coupled- 
cluster  wave-function  is  used  as  the  ground  state  wave-function,  the  one-hole  and 
one-particle  sector  decouple  with  respect  to  the  eigenvalues.  The  self-energy  is  then 
shown  to  be  equivalent  to  the  correlation  potential,  derived  from  the  partitioned 
equation-of-motion  coupled-cluster  method. 

In  Chapter  4 the  infinite-order  correlation  potential,  derived  from  the  Fock  space 
coupled-cluster  approach,  is  approximated  by  a potential  correct  through  second 
order  in  perturbation.  Since  the  one-hole  and  the  one-particle  sectors  are  decoupled 
in  the  Fock  space  coupled-cluster  method,  the  orbitals,  obtained  from  the  CIP  model, 
are  rotations  within  the  occupied-occupied  space  and  rotations  within  the  virtual- 
virtual  space  and  are  therefore  not  an  improvement  over  Hartree-Fock  orbitals.  To 
improve  the  orbitals,  an  additional  condition  has  to  be  applied.  The  Brillouin- 
Brueckner  condition  is  used  to  assure,  that  the  single  determinant,  that  is  formed 
by  the  orbitals,  has  maximum  overlap  with  the  exact  wave-function.  The  Brueckner 
condition  is  ensured  through  second  order  in  perturbation. 

Two  different  partitioning  schemes  of  the  Hamiltonian  were  used,  leading  to 
two  different  second-order  Hamiltonians.  Results  for  ionization  potentials,  electron 
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affinities,  and  dipole  moments  are  given  for  some  sample  molecules.  The  dipole 
moments  show  significant  improvement  over  the  Hartree-Fock  values,  if  the  second- 
order  Brueckner  determinant  is  used. 

The  second-order  IP’s  and  EA’s  show  mostly  an  improvement  over  Koopmans’ 
values.  However,  it  became  apparent  that  intruder  state  problems  are  severe;  in 
certain  cases  intruder  states  from  the  shake-up  eigenspectrum  even  occur  among 
the  principal  IP’s.  Certainly,  the  quality  of  the  IP’s  and  EA’s  is  influenced  by 
intruder  states.  In  Section  3.4.2  it  was  shown,  how  to  obtain  Fock  space  amplitudes 
through  EOM-CC  amplitudes.  In  order  to  ensure,  that  intruder  states  do  not  occur 
in  the  principal  IP  and  EA  spectra  of  the  second-order  effective  Hamiltonian,  the 
EOM-CC  problem  could  be  expressed  up  to  second  order  in  perturbation,  and  the 
second-order  Fock  space  amplitudes  could  be  calculated  through  the  second-order 
EOM-CC  amplitudes,  as  has  been  shown  for  the  infinite-order  expressions. 
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